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Synopsis 
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We present h-p type Spectral Element Methods for Elliptic Boundary Value Problems 
on Polygonal Domains using Parallel Computers. We first discuss the Poisson 
equation on a polygon. 

To resolve the singularity which arises at the corners we use a geometric mesh as 
has been proposed by Babuska and Guo. With this mesh we seek a solution which 
minimizes a weighted norm of the residuals in the partial differential equation and 
a fractional Sobolev norm of the residuals in the boundary conditions and enforce con- 
tinuity by adding a term which measures the jump in the function and its derivatives 
at inter-element boundaries, in an appropriate fractional Sobolev norm, to the func- 
tional being minimized. Since the second derivatives of the actual solution are not 
square integrable in a neighborhood of the corners we have to multiply the residuals in 
the partial differential equation by an appropriate power of r*,, where measures the 
distance between the point P and the vertex Ak in a sectoral neighborhood of each of 
these vertices. In a neighborhood of the corner Ak we switch to new variables {Tk,0k) 
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where r* = Inr^ and {rkyOk) are polar coordinates with origin at Ak- In doing so the 
geometrical mesh is reduced to a quasi-uniform mesh in a sectoral neighborhood of 
the corners and so Sobolev’s embedding theorems and the trace theorems for Sobolev 
spaces apply for functions defined on mesh elements in these new variables with a 
Uniform Constant. We then derive differentiability estimates with respect to these 
new variables and a stability estimate for the functional we minimize. 

To solve the minimization problem we have defined, we need to solve the normal 
equations for the least-squares problem corresponding to collocating the partial differ- 
ential equation and the boundary conditions at an over-determined set of collocation 
points and enforcing continuity of the function and its derivatives at the collocation 
points at inter-element boundaries, suitably weighted. However, we do not need to 
compute and store any matrices, like the mass and stiffness matrices, to compute the 
residual in the normal equations. 

We then use the stability estimates to obtain parallel preconditioners and error 
estimates for the solution of the minimization problem. We can precondition the normal 
equations by using this preconditioner which is of block diagonal form and which allows 
the solutions for different elements to decouple completely. Moreover this is nearly 
optimal as the condition number of the preconditioned system is polylogarithmic in 
N, which is proportional to the number of processors and the number of degrees of 
freedom in each variable on each element. Finally we show that the error we commit 
is exponentially small in N. 

For the purely Dirichlet problem our spectral element functions are non- conforming 
and hence there are no common boundary values to solve for. This no longer holds 
for problems with Mixed boundary conditions. Here our spectral element functions 
are essentially non-conforming except that they are continuous at the vertices of the 
elements on which they are defined. Hence the set of common boundary values are 
the values of the function at the vertices of the elements and so their cardinality, 
being proportional to the number of elements, is small as compared to Finite Element 



Methods, where the common boundary values are the values of the functions along the 
edges of their elements. Thus the method proposed is a Vertex-Based Method. 

We solve the normal equations by the preconditioned conjugate gradient method. We 
first solve a much smaller system of equations corresponding to the Schur complement 
of the sub- vector of common boundary values for which we need to be able to compute 
a preconditioner. The Schur complement matrix can be computed accurately since 
it’s dimension is small. This method turns out to be computationally more efficient 
than Finite Element Methods where complex techniques have to be used to obtain a 
preconditioner for the Schur complement matrix. Moreover since for the problem with 
Dirichlet boundary conditions the spectral element functions are non-conforming, there 
is no set of common boundary values and so the method is even more efficient. 

We should mention that once we have obtained our approximate solution consisting 
of non-conforming spectral element functions we can make a correction to it so that 
the corrected solution is conforming and is an exponentially accurate approximation 
to the actual solution in the norm over the whole domain. 

All these results are valid for elliptic problems with mixed boundary conditions on 
domains with curvilinear boundaries which satisfy the usual conditions. 
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Chapter 1 


Introduction 

This dissertation is on h-p Spectral Element Methods for Elliptic Problems 
on Non-smooth Domains using Parallel Computers. We consider here the two 
dimensional problem. 

1.1 Review of Existing Work 

Spectral methods are one of the three major techniques for the numerical solution of 
Partial Differential Equations (PDE), namely Finite Difference Methods (FDM), Fi- 
nite Elements Methods (FEM) and Spectral Methods (SM). Among these, Spectral 
methods have emerged as among the most accurate way of solving PDE and have 
become increasingly popular for problems where high accuracy is desired for complex 
problems like numerical simulations of turbulent flows, numerical weather prediction, 
ocean dynamics etc. If properly constructed these methods converge to the solution 
with exponential accuracy and can be used to obtain solutions where other numerical 
techniques fail. The very high accuracy of Spectral methods allows one to treat prob- 
lems which would require an enormous number of mesh points by finite differences, 
with much fewer degrees of freedom. However, a poorly designed Spectral method may 
perform much worse than simpler FDM. 

In 1970s Gottlieb, Orszag and others initiated the work on problems in fluid dynam- 
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ics using Spectral methods and presented the formulation of modern Spectral methods 
in the monograph [22]. The book of Canute, Hussaini, Quarteroni and Zang [14] fo- 
cuses on fluid dynamics algorithms and includes practical as well theoretical aspects 
of Spectral methods. Among other references are the book by Voigt, Gottlieb and 
Hussaini [45] and the latest one in this series is by Karniadakis and Sherwin [28]. 

Spectral methods are based on global interpolants, in contrast to FDM and FEM. 
The basic principle of this numerical technique is to represent the dependent variable 
in a finite (truncated) series of known infinitely difiFerentiable global functions, for 
instance, 

N 

U (x) ~ UN (x) = <^n<Pn, 

71=0 

where (x) are the approximation functions (also known as the trial or expansion 
functions) and which are used as the basis functions for a truncated series expansion 
of the solution. The series is then substituted into the differential (or integral) equa- 
tion and upon the minimization of the residual function the unknown coefficients are 
computed. 

Most commonly used approximation functions are Trigonometric polynomials (for 
periodic boundary conditions), Chebyshev polynomials and Legendre polynomials. The 
three most commonly used methods are Collocation^, Galerkin and Spectral Tau, which 
differ in their choice of test (or weight) functions. The Collocation approach is the 
simplest and most widely used among these methods. The most effective choice for the 
grid points are those which corresponds to quadrature formula of maximum precision, 
like Legendre-Gauss-Lohatto points. 
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1.2 Spectral Methods on Non-smooth Domains 

Unlike FEM and FDM, the order of the convergence of Spectral methods is not fixed 
and it is related to the maximum regularity (smoothness) of the solution. The high 
accuracy of Spectral methods is easily lost if the solution has finite regularity (i.e. in 
the presence of discontinuities) or the domain is irregular. 

The latest book on Spectral/fip Element Methods for CFD by Karniadakis 
and Sherwin [28] summarizes the recent research in the subject. Here we present a brief 
summary of the Section Non-Smooth Domains in the Chapter Helmholtz Equation. 

Exponential convergence is obtained with spectral/hp element methods if the so- 
lution is smooth, possessing a high degree of regularity. However, there are a number 
of cases for which the solution of a Poisson’s equation may be singular, like solutions 
in a non-smooth domains, smooth domains but with a discontinuity in the boundary 
conditions or in the specified data (eg. forcing). Here, as in [28], we also assume that 
all the data, as well as the boundary conditions, are smooth and that singularities 
are due to corners in the domain. First derivatives are unbounded when the angle is 
reflexive, and second derivative are unbounded when the angle is acute or obtuse. In 
this case not only the fast convergence of spectral/ ftp discretizations may by lost, but 
also the numerical solution obtained (with any standard method) may be erroneous. 

There are four methods discussed in that section, that allow us to recover, if not ex- 
ponential, at least very fast convergence for most elliptic problems. These are described 
in the next few subsections. 

1.2.1 Gradual h Refinement or h/hp FEM 

This approach requires good discretization and quasi-uniform meshing. A geometric 
mesh has been found to be effective with a ratio of 0.15 [11]. The method gives 
exponential convergence. 

The global matrix A can be split into components containing boundary and interior 
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contributions, that is, 


An Aib 
Abb 


where Abb represents the components of A resulting from boundary-boundary mode in- 
teractions, Aib represents the components resulting from coupling between boundary- 
interior modes and An represents the components resulting from interior-interior mode 
interactions. 

In the global system the global boundary-boundary, Abb, matrix is sparse and 
may be reordered to reduce the bandwidth. The elements of Abb, corresponding to 
the common boundary values, are the values along the edges of the elements. 

To solve the system Au = h the block L- U factorization of A 


I 0 



0 


I Ajl Aib 

1 

1 


0 

S 


0 I 


where the Schur complement matrix S is defined as S' = Abb — AJ^AJ^Aib, is 
used. Consequently, solving Au = h reduces to solving Sub = Lb, where Lb = 
Hb — AjgAjlhi. Once - ub is known we can determine ui from the equation ui = 
■^ 7 / ~ -^ibUb- 


Preconditioning 

To solve the algebraic system of equations arising from the discretization of symmetric 
elliptic BVPs via spectral/hp methods we need very efficient and effective precondi- 
tioners for iterative solvers for large scale computations in parallel environments. 

The iterative solution techniques for the linear systems generated from the h- version 
or p- version of the FEM have been widely studied in the past decades. Among them 
most successful are the preconditioned conjugate gradient methods. The efficiency of 
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contributions, that is, 



Aji Ajb 
(^/b)^ Abb 


where Abb represents the components of A resulting from boundary-boundary mode in- 
teractions, Aib represents the components resulting from coupling between boundary- 
interior modes and An represents the components resulting from interior-interior mode 
interactions. 

In the global system the global boundary-boundary, Abbi matrix is sparse and 
may be reordered to reduce the bandwidth. The elements of Abb, corresponding to 
the common boundary values, are the values along the edges of the elements. 

To solve the system Au = h the block L- U factorization of A 





I 

o 

1 



0 s 



I Ajj Aib 
0 I 


where the Schur complement matrix S is defined as S' = Abb — ^b-^ii-^ib, is 
used. Consequently, solving Au = h reduces to solving Sub = Lb, where Jib = 
Hb — nhn Once wB is known we can determine uj from the equation uj — 
Ajj hi — Ajj AibUb- 


Preconditioning 

To solve the algebraic system of equations arising from the discretization of symmetric 
elliptic BVPs via spectral//ip methods we need very efficient and effective precondi- 
tioners for iterative solvers for large scale computations in parallel environments. 

The iterative solution techniques for the linear systems generated from the h- version 
or p- version of the FEM have been widely studied in the past decades. Among them 
most successful are the preconditioned conjugate gradient methods. The efficiency of 
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these solution techniques relies largely on the condition number of the preconditioned 
system. For the h- version, Bramble et al. [12] proposed a class of preconditioners 
based on a representation of the H^^— norm on the interfaces of subdomains, which 
generally results in a condition number of order O (1 + In , where H and h are 
the diameters of the subdomains and the elements respectively. On the other hand, 
Babuska et al. [5] proposed a preconditioner for the p- version of the FEM, based 
on the orthogonalization procedure. The resulting condition number is of the order 
O (1 + Inp)^ with p being the polynomial degree used. 

The preconditioning for the h-p version has been considered only recently. In the 
paper of Guo and Cao [26], a preconditioner for the h-p version of the FEM in two 
dimensions, with non-uniform mesh and non-uniform distribution of element degrees 
has been proposed which is also applicable to h-p version associated with the geomet- 
ric mesh. The condition number of the preconditioned linear system is of the order 
maxj ^1 -f- In , where Hi is the diameter of the subdomain fij, hi and pi are the 
diameter of elements and the maximum polynomial degree used in flj. This fully covers 
the h- version when Pi = 1 and p- version when hi = Hi as special cases. 


1.2.2 Conformal or Auxiliary Mapping 

We will consider three different situations for individual equations, namely Laplace, 
Poisson and Helmholtz equation. 


Laplace Equation 

For the Laplace equation Au = 0, with homogeneous boundary conditions the solution, 
in the neighborhood of the corner, can be expressed as 

OO 

u{r,e) = Y^ak(j)k{r,0) 

k=0 


5 
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where the coefficients are determined by the boundary conditions and 


<f>k {r,0) 


sin {k^O) , if not an integer 

[in r sin (k^O) + 6 cos {k~Gy\ , if A: ^ is an integer 


Here oj is the sectoral angle. Assuming that the logarithmic term doesn’t contribute to 
the solution, the mapping z = which is conformal at all points except the origin 
makes the transformed solution analytic in terms of the new variables. Thus 


u (r, 9) = sin (^k—9^ 

k=o ^ 


OO 

u (p, •!/;) = ^ akp^ sin [k^)) . 

fc =0 


This method was first implemented in [10] for hp finite element discretizations. 


Poisson Equation 

For the Poisson equation Au = f{x), the situation is more complicated because the 
solution may not be analytic after the mapping. Typically we decompose the solution 
into a homogeneous part u^, which has the singularity, and a particular part which 
depends on the forcing. Complications arise due to the particular solution since even 
if it is smooth in the original domain it may be singular after the transformation. 

In order to enhance the convergence, we separate the two contributions so that we have 
an analytic contribution in the z plane and an analytic contribution in the ^ plane. 


Helmholtz Equation 

For the Helmholtz equation Att — Xu = f (x) , A > 0, the conformal mapping is an 
effective way of enhancing convergence although exponential convergence cannot be fully 
restored. The auxiliary mapping z,= converts the Helmholtz equation to 

Au - A ^ 
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In terms of the original variables the solution around the corner is 

CX) 

u (r, 0) = ^2 ^\/Ar^ sin 

for not an integer, where Im {z) is the modified Bessel function of first kind. After 
application of the mapping the solution has the form 

oo / oo 

w (p, '0) = X) 

fc=o Vi=o 

with leading singular term of order . Therefore, the estimated convergence rate 
is O , which is algebraic but in practice adequately fast. 

1.2.3 Singular Basis 

An alternative approach to using the auxiliary mapping is to use a set of supplementary 
basis function which have the leading behavior of the singularity in conjunction with 
the smooth basis {x) ■ For the Helmholtz equation the leading order singular terms 
are 

r^, (->1/2), (->l),..., 

TT \7r / 

which can be included into the expansion basis. However, we can do even better by 
supplementing the standard basis in the mapped domain. The transformed solution is 
then 

OO OO OO 

u (p, ■0) = ^ ttkhi (y^P ^ ) sin k-tp 

k=l k=l 1=0 

and thus the leading singularities are weaker, i.e. 

pl+2^, p2+2H 0>l/2), 

It has been shown in [39] that with one or two terms included in the transformed 
domain, a ery fast convergence is obtained. To achieve the exponential accuracy 
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we need to include higher order terms but in the limiting case the system becomes 
ill-conditioned. 

1.2.4 Steklov Formulation 

A more recent method that treats singular solutions of both scalar and vector elliptic 
problems in the neighborhood of corners has been suggested by Szabo and Babuska. As 
we have already seen such singular solutions are characterized by the form u = r^f (9) 
close to the corner. We consider the given domain and choose the radius of the domains 
R sufficiently small so that the general form of the solution is still valid and thus 

1 

— = j3r^-^f {$) = (p/R) u on Ts 

where Fa denotes the circular arc of radius R. We consider the Laplacian equation with 
Robin boundary conditions on the segments, which enclose the corner, Fi and F2, that 
is, 

au {x) + =0, 6^0. 

on 

The variational statement of the Laplace equation with the above boundary conditions 
is as follows: 

Find real numbers /? and u € H^ {Or) so that 

2 

a {u, v) + ^ Mi {u, v) = PMz {u, v) Vu G H^ 

2=1 

where a{u,v) = (Vu, Vu) is the standard weak Laplacian operator and we have 

Mi (it, v) = ^ 

We now extend the variational formulation of the problem to a domain that does not 
include the singular vertex. To this end we consider a cut in the form of circular arc 
R*, denoted by F4, that surrounds the corner. The new domain is smaller than before 
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and doesn’t contain any singularities. The modified Steklov problem has the form: 
Find real numbers and u € H^ (Or*) so that 


a {u, v) + Mi {u, v) = fi 


Z=1 


^ Mi {u, v) 


i=2 


VveH^(nR^), 


where M 4 (u, v) is defined as the rest of Mi {u, v) , (i < 3) evaluated on r 4 . Since the 
domain Qr» does not contain any singularities, standard spectral//ip discretization 
can be used to obtain spectral convergence. This method is set in the finite element 
formulation. 


1.3 A Review of the Thesis 

Current formulations of spectral methods to solve elliptic problems in non smooth 
domains allow us to recover only algebraic convergence [14, 28]. One method, which 
yields relatively fast convergence, makes use of a conformal mapping of the form z = 
to smooth out the singularity that occurs at the corner and is referred to as the method 
of Auxiliary Mapping. However “even though the conformal mapping is an effective way 
of enhancing convergence, exponential convergence cannot be fully recovered" [28]. 

A method for obtaining a numerical solution to exponential accuracy for elliptic 
problems with analytic coefficients posed on curvilinear polygons whose boundary is 
piecewise analytic with mixed Neumann and Dirichlet boundary conditions, was first 
proposed by Babuska and Guo [7, 8] within the framework of the FEM. They were able 
to resolve the singularities which arise at the corners by using a geometrical mesh. 

We also use a geometrical mesh to solve the same class of problems to exponential 
accuracy using h-p Spectral element methods but with an important difference. With 
this mesh we seek a solution which minimizes a weighted norm of the residuals in 
the partial differential equation and a fractional Sobolev norm of the residuals in the 
boundary conditions and enforce continuity by adding a term which measures the jump 
in the function and its derivatives at inter-element boundaries, in an appropriate frac- 



tional Sobolev norm, to the functional being minimized. Since the second derivatives of 
the actual solution are not square integrable in a neighborhood of the corners we have 
to multiply the residuals in the partial differential equation by an appropriate power of 
ffc, where rk measures the distance between the point P and the vertex Ak in a sectoral 
neighborhood of each of these vertices. In each of these sectoral neighborhoods we use 
a local coordinate system {rk,6k) where = Inrfc and {rk,9k) are polar coordinates 
with origin at Ak- We then derive differentiability estimates with respect to these new 
variables and a stability estimate for the functional we minimize. All these quantities 
are computed in different coordinate systems consisting of local coordinate systems in 
a sectoral neighborhood of each of the vertices and a global one outside these neighbor- 
hoods. Thus we use the auxiliary map 2 : = log^ to solve the problem with exponential 
accuracy. 

We then use the stability estimate to obtain parallel preconditioners and error esti- 
mates for the solution of the minimization problem which are nearly optimal as the con- 
dition number of the preconditioned system is poly logarithmic in N, i.e. O (1 -t- In 
where N is the number of processors and the number of degrees of freedom in each 
variable on each element. Moreover if the data is analytic then the error is exponentially 
small in N. 

The method we propose is essentially a least-squares method. To compute the resid- 
uals in the normal equations, however, we do not need to compute any mass and stiff- 
ness matrices nor do we need to filter the coefficients of the differential and boundary 
operators or the data. We solve the resulting normal equations by the preconditioned 
conjugate gradient method. 

To solve elliptic boundary value problems with mixed Neumann and Dirichlet 
boundary conditions the spectral element functions we use are not fully non-conforming; 
instead the spectral element functions are restricted so that their values at the com- 
mon vertices of the elements on which they are defined are the same. We divide the 
vector composed of the values of the spectral element functions at the Legendre-Gauss- 
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Lobatto points on each element into two sub vectors - one consisting of the values of the 
spectral element functions at the vertices of all the elements and the other consisting 
of the remaining values. The dimension of the first sub vector, which can be thought 
of as a vector of common boundary values is very small compared to FEM. 

The normal equations thus obtained are solved by the preconditioned conjugate 
gradient method. We first solve a much smaller system of equations corresponding 
to the Schur complement of the sub vector of common boundary values. The Schur 
complement matrix can be computed accurately since it’s dimension is small. Moreover, 
the error in the numerical solution is exponentially small in N, which is proportional 
to the number of elements and the degree of the spectral element representations in 
each variable on each element. 

We finally show that the h-p Spectral element method applies to elliptic problems on 
curvilinear polygons with mixed Neumann and Dirichlet boundary conditions provided 
the Babuska-Brezzi inf-sup conditions are satisfied. 

We now briefly describe the contents of this dissertation. In Chapter 1 we provide 
a review of existing work. In Chapters 2 and 3 we restrict ourselves to examining the 
Poisson equation with Dirichlet boundary conditions on a polygon. In Chapter 2 we 
obtain differentiability estimates in modified polar coordinates and prove the stability 
theorem 2.3 on which our method is based. Since the statement of this theorem may 
appear complicated we try and provide motivation for it by stating the stability theorem 
2.1 for a simpler case. In Chapter 3 we provide computational techniques and error 
estimates for Dirichlet problems. 

In Chapter 4 we examine how to solve the Poisson equation with mixed Dirichlet and 
Neumann boundary conditions. In both Chapters 3 and 4 we also provide numerical 
results. 

Finally in Chapter 5 we generalize all our results to elliptic problems with analytic 
coefficients, posed on curvilinear polygons with piecewise analytic boundaries, which 
satisfy the Babuska-Brezzi inf-sup conditions. 



Chapter 2 


Differentiability and Stability 
Estimates for Dirichlet Problems 

2.1 Introduction 

To overcome the singularities that arise in a neighborhood of the corners we use a 
geometrical mesh, as has been proposed by Babuska and Guo. The geometrical mesh 
becomes geometrically fine in a neighborhood of each of the corners. In a neighborhood 
of the corner Ak we switch to new variables {Tk,6k) where Tk = Inrfe and (rk,9k) are 
polar coordinates with origin at Ak- In doing so the geometrical mesh is reduced 
to a quasi-uniform mesh in a sectoral neighborhood of the corners and so Sobolev’s 
embedding theorems and the trace theorems for Sobolev spaces apply for functions 
defined on mesh elements in these new variables with a Uniform Constant. These 
new variables, which we shall refer to as modified polar coordinates, were first used 
by Kondratiev in his seminal paper [30]. Away from these sectoral neighborhoods of 
the corners we retain {x, y) variables for our coordinate system. Thus we also use the 
auxiliary map z = log^ to remove the singularity at the origin and this enables us to 
obtain the solution with spectral accuracy. 

By subtracting an analytic function from the solution if necessary, we may assume 
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that the Dirichlet data vanishes at the corners. We seek an approximate solution which 
vanishes at the corner-most elements and is a sum of tensor product of polynomials of 
degree N in and Ok in the remaining elements of the sectoral neighborhood of the 
corners. The remaining quadrilateral elements are mapped to the unit square S and 
the approximate solution is represented as a sum of tensor products of polynomials 
of degree iV in ^ and 77 , the transformed variables. If Neumann boundary conditions 
are imposed on both the sides which meet at the corner, the approximate solution at 
corner most elements is represented by a constant, instead of zero. 

We now seek a solution as in [17] which minimizes the sum of the squares of a 
weighted norm of the residuals in the partial differential equation and the sum of 
the squares of the residuals in the Dirichlet boundary conditions in an appropriate 
Sobolev norm and enforce continuity by adding a term which measure the sum of 
the squares of the jump in the function in the norm and the squares of the jump 
in its derivatives across inter-element boundaries in an appropriate fractional Sobolev 
norm to the functional being minimized. Since the residuals in the partial differential 
equation blow up in a neighborhood of the corners, we have to multiply these residuals 
by an appropriate power of rfe, where r*, measures the distance between the point P 
and Ak- All these computations are done using modified polar coordinates in a sectoral 
neighborhood of the corners and a global coordinate system elsewhere. 

For the Dirichlet problem we use Spectral element functions which are non- conforming. 
To solve the minimization problem we have defined, we need to solve the normal equa- 
tions for the least-squares problem corresponding to collocating the partial differential 
equation and boundary conditions at an over- determined set of collocation points and 
enforcing continuity of the function and its derivatives at the collocation points at 
inter-element boundaries, suitably weighted. However we do not need to compute and 
store any matrices, like the mass and stiffness matrices, to compute the residual in the 
normal equations [18]. 

We can precondition the normal equations by using a preconditioner which is of 
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that the Dirichlet data vanishes at the corners. We seek an approximate solution which 
vanishes at the corner-most elements and is a sum of tensor product of polynomials of 
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corners. The remaining quadrilateral elements are mapped to the unit square S and 
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equation blow up in a neighborhood of the corners, we have to multiply these residuals 
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and Afe. All these computations are done using modified polar coordinates in a sectoral 
neighborhood of the corners and a global coordinate system elsewhere. 

For the Dirichlet problem we use Spectral element functions which are non-conforming. 
To solve the minimization problem we have defined, we need to solve the normal equa- 
tions for the least-squares problem corresponding to collocating the partial differential 
equation and boundary conditions at an over- determined set of collocation points and 
enforcing continuity of the function and its derivatives at the collocation points at 
inter-element boundaries, suitably weighted. However we do not need to compute and 
store any matrices, like the mass and stiffness matrices, to compute the residual in the 
normal equations [18]. 

We can precondition the normal equations by using a preconditioner which is of 
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block diagonal form and which allows the solutions for different elements to decouple 
completely. Moreover this is nearly optimal as the condition number of the precondi- 
tioned system is polylogarithmic in N, which is proportional to the number of proces- 
sors and the number of degrees of freedom in each variable on each element. Finally 
we show that the error we commit is exponentially small in N. 

2.2 Function Spaces and A’Priori Estimates 

Let be a polygon with vertices Ai, A 2 , . . . ,Ap and corresponding sides Fi, r 25 . - . , Fp 
where Fj joins the points Ai_i and Aj. In addition let the angle subtended at Aj be Uj. 



Figure 2.1; Polygonal domain 


In this chapter we shall examine the solution of the problem 

(2.1a) Au = / for {x, y) € Q, 

with Dirichlet boundary conditions 

(2.1b) u = Qj for (a:,y) 6Fj, or 

u = for {x, y) e F^°^ = 90. 

Let Z denote the point Z = {x,y) . 

We now need to review a set of a’priori estimates proved in {?]. Let (0) denote the 



2.2. FUNCTION SPACES AND A’PRIORI ESTIMATES 


15 


completion of the space of infinitely differentiable functions with respect to the norm 


|u 


2 _ 
m,n 



dxdy. 


Let Pi denote the Euclidean distance between Ai and Z, i.e. pi = \Z — Ai\-, we then 
define rj = min(l,pi). We shall let ^ denote the multi index § = ,Pp)- 

Further, we define 


(Z) = (Z) . 


By (Q) we denote the completion of infinitely differentiable functions with respect 
to the norm 

m 

k=l^\a\=k 

Let |.|jg gpace of functions <f>j such that there exists / € (f7) 

so that /Ifj- = and define 


ll<^. 




m-l/2,i-l/2,„ , = 


e 


(r,) 


inf 


11/1Ih; 




Let 


4 (f2) = {u {Z) |w G (0) , m > i I 


and 


QSi, (fi) = 


U(Z) ue (Q) , ||1£|“«| «s+»-,|li.(„) < Cd’-‘ 
for |q;| = A: = Z, 1 4 - 1 , . . . , d > 1 , C independent of k 


Let Q C K.^ be an open set with a piecewise analytic boundary dQ and let 7 be part 
or whole of the boundary dQ. Finally let ( 7 ) , 0 < Z < 2 , be the space of all 
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functions (p for which there exists / E 58^ (Q) such that / = <^ on 7. 

We now cite the important regularity theorem 2.1 of [7]. Let / € E 

553/2 ^p[o]^ ^ 

P = {^ 1 , P 2 , , Pp) , 0 < Pi < 1, Pi > 1 — 'njiOi. Then the problem (2.1a - 

2.1b) has a unique solution in (f2) and u E 55^ (0) . 

Now in Section 4 of [8] it has been shown that when is analytic on every closed 
arc Ti and is continuous on Ft®] then € 05^^^ (pM) . Further if / is analytic then 
it belongs to 58^. 

Next as in [8] we introduce the space : 

I ue Hf (Q) \\D-u{Z)\ < Od'^kl (Z))-' , 

I |q:| = A: = 1, 2, . . . , C > 1, d > 1 independent of k 

The relationship between and 58^ is given by Theorem 2.2 of [8] which we state as 
follows: 



^ m C €% 

Finally we need one last result from [8], viz. Lemma 2.1 which is stated below: 

Let u E (fi) . Then u is continuous on Q, and 

ll^llc(n) — ^ ll'^llif|’^(a) • 

Since we are assuming that the data, gi,... ,gp are analytic and compatible at the 
vertices the values of u at the vertices AijA^,... ,Ap are well defined. Thus if we 
subtract from u an analytic function which assumes these values at the vertices then 
the difference would satisfy (2.1a - 2.1b) with a modified set of analytic data and the 
Dirichlet boundary data would assume the value zero at Ai, A2^ ■ . ■ , Ap. Hence without 
loss of generality we may assume gj-i (Aj) = gj (Ay) = 0 for j > 1 and gp (Ai) = 0. 
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We now define one last norm which will be needed in the sequel. Let 


xl4, Ann 


2 rVii r 




Dr-D^^U 


dTjdOj. 


pipu n\nti 2 

(2.2) E / . / (v»?0 

J-oo V ^ W 


We need to obtain an asymptotic estimate on as /x — > 0. 

It is easy to estimate the terms in the right hand side of the above definition when 
|a| > 2. For 

f^ipi AnfjL 

/‘Vi rlniii . . 2 

< dTjdOj 

r'i’i rM 

2<|a|<m'^’^ •'0 


^2(1/5,) ^ rjdr. 


jdOj 


Here we have used the fact that m G (S^) and Theorem 1.1 of [7] to obtain the above 
result. Next we bound the terms when |a| = 1. Since u € (Q) we have 


(2.3) 


\D^u {Z)\ < Cd {Z))-^ when la| = 1. 


Moreover we have 


/ {^Tj + ^e,) dTjdOj = + ^y) dxdy. 


Here 


Sj = {(a:,y) |0 < rj < fx, < Oj < •0^} . 


Hence by the above relation 


(2.4) J j (it| + u^) dxdy < 2C^d^ J j Tj "^^Wjdrjddj < {Rdf 
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Finally we have to estimate 



u ( tj , 6j)f dTjddj. 


Since u vanishes at Aj 



iv, Oj) drj. 


Hence 




[ Up (p, Oj) dp 

Jo 


Here p = e'' and rj = . Since u € (0) we obtain 




And integrating the above with respect to tj and 9j gives 


(2.5) 


fxpi pin 11 

/ / \u{Tj,6j)f dTjdOj < [Kd^ 

J ib^j J —OO 


flpj J —OO 


Combining (2.2), (2.4) and (2.5) we get the required estimate 


(2,6) 


ipi pin II 


^ /; / 


ja]<m 


Ut^ 10 J 2 


dTjddj < (Cd”‘-2(m-2)!)^ 


2.3 Stability Estimates 

We first need to divide Q into subdomains. Thus we divide Q into p subdomains 
S^,S'^,. . . ,S^, where 5* denotes a domain which contains the vertex Ai and no other, 
and on each 5* we define a geometric mesh as has been done in [8]. 
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Figure 2.2: Mesh over the whole domain 
Ufcrri Then © satisfies the following conditions: 

1- are curvilinear quadrilaterals or triangles and the intersection of any two j 
is one common vertex or one entire side or is empty. 

2. Let and be the maximal and minimal length of the sides of We shall 
assume there is a constant independent of i, j, k and of the partition such that 

(2.7) < A. 

3. Let M = {M^j, 1 <i < hj, 1 < j < 1 < fc < p} in which M^j is a one-to- 

one mapping of the closed standard master square S = [0, 1] x [0, 1] {respectively 
standard master triangle T = {(^, p)|0<?7<l — 0<^<1} onto Let 

P^j i and denote the vertices and sides of then [MiJ) ^ 

(AIij)~^ (jlj^i) denote the vertices and sides of S (respectively T), 1 < Z < 4 
(respectively 1 < Z < 3). Moreover if and map the closed standard 

I 

square S onto elements and ^ with common side 7 = P1P2, then for any 
P e 7. dist (P), (P,)) = dist (P), (P,)), 
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1 < i < 2. We will assume M/',- can be written in the form 

(2.8) X = (f, >;)€ S (respectively T) 

V = 

with Xi j and being analytic functions on S (respectively T). Further we 
assume that for jal < 2 

(2.9a) \D‘^x\,\D^y\ < and 

(2.9b) 4 . <C,(hyy 

for all i,j and k with constants C,CiC 2 independent of i,j and k and being 
the Jacobian of the mapping . 

Let fj, = {fxi, . . . , fip) with 0 < < 1. Then 6^ is called a geometrical mesh with 

ratios /lx = (/xi, . . . , fip) when in addition the following condition is fulfilled. 

4. Let ,• e © and dix denote the distance between ,■ and A*.. Then and h'-x 
satisfy 

(2.10a) Cl < 4 . < C 2 ,l<j<N, 

(2.10b) Czp < d^j < C^p, N <j <Jk, l<i< Ik,j, 

(2.10c) dl, = 0, 1 < i < Ik,i, l<k<p, 

(2.10d) K,dl^ < hlj < hlj < K2dl^ 

for 1 < 7 < Jfc, 1 < i < Ik,ji 1 < ^ < P) where C; for 1 < 1 < 4 and Ki for 
1 < Z < 2 are constants independent of i, j and k. Moreover J}. = N + 0(1). 

We now put some restrictions on 6. Let (rky^k) denote polar coordinates with center 
at Ak- Let Tk = Inr*;. We choose p so that the sector S* with sides Tk and F^+i, center 


1 ^ ^ Ik,jj 

1 < k < p, 
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at Ak and radius p satisfies 


sjs U 


n}’ -es* 


may be represented as 


{2.11a) 


Sp = {( 3 ^> y) |0 <rk<P, ipf <ek<ip^}. 


Let be an increasing sequence of points such that 

ip!^. Let Aipf = 7/)*+]^ — ip^. We choose these points so that 


k _ 


(2.11b) 


max ^max < A m^in ^min j 


for some constant A. Let 


(2.11c) 

(2.11d) 


(Tj = 0, and 

aj = p (//jfc)^'^^'"^ for 2 < j < iV + 1. 


Finally we define 

( 2 . 12 ) 


rj’^ =lna^foTl<j <N + 1. 


Let 


(2.13) 


y) ki <^k< cTj+i, < Ok < i^i+i } , 


for 1 < z < /jt, 1 < ^’ < N. 

We assume there exists a number u such that 


i’l and = 


(2.14) 


^i,N+i = y) |p < ■rjfc < V’,* < Ok < ipi+i } , 
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for 1 < z < Ik- 

In other words hj is independent of j for j < N +1. We shall let O’" denote 

(2.15) O’" = |1 < z < 4, 1 < j < W} , 

for 1 < A: < p. 

Let 


(2.16) 0^+' = {Of,,. |1 < z < N + l<j<Jk,l<k<p}. 

We shall relabel the elements of 0^+^ and write 


(2.17) OP+^ = {Qf+^|l</<L}, 

where L denotes the cardinality of We shall let Q’" denote the sector with vertex 
at Ak given by 


(2.18a) 

(2.18b) 


O’" = {(a;, y) |0 < < p, tpf < dk < 'ip^} , and 

= q\ j|Jo''|. 


Notice that all the sets Q.’" are open sets. 


Henceforth to keep our notation simple we will assume that are quadrilaterals 
ioT 1 < k < p, N + I < j < Jk, I < i < h- Moreover we assume Ikj < I for all k and 
j. Here J is a small integer, and this fact plays a fundamental role in allowing us to 
use non-conforming spectral elements to solve the Dirichlet problem. Further let 
denote the sides of the quadrilateral 1 < A < 4. Then we assume 


(2.19) 



y = 


1 = 1,3 
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( 2 . 20 ) 



y = 


0 < 7? < 1, i = 2,4 


where 4>Xi,v V’fj.z analytic functions for all i,j, k and 1. We wish to obtain a stability 
estimate for the function u, given by a non-conforming finite dimensional representation 
ufj on each domain for the entire polygonal domain Cl. Now, as stated in the 
introduction we partition the open set Cl into p open sets S^,S^,... ,S^ such that 
each S'* contains only the singularity at the vertex Aj. Let S’^ be one of these open 
sets. Then S^ = (J IJ T*' where Cl^ is the open sector with center at A^ and 
radius p, is the circular arc which bounds fl*, and T^ is the open set defined as 

r'' = S'=\(f2^U^p) • 



The domain S* is as shown in Fig. 2.3. Two of its sides are the straight lines 
Ck+iOdS^ and F^pl^S*^. The remaining side dS^ consists of piecewise analytic arcs. 
The subscript c in dS^ denotes curvilinear. S* is partitioned by a set of arcs { 7 /}^ into 
subdomains. 


= {(^> y) ki <n< cTj+i, <0k< ^'i+i } 


Now 
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for 1 < i < /fe, I < j < N. Let Tk = Inrfc. Define 


( 2 . 21 ) 


= {in, h) \nj <Tt< 7)j+i, V’* < < lAj+i } 


for 1 < i < Jjt, 1 < j < N. Let w be a smooth function. Then 


i^xx 


I f d d(jj \ -2t / 


Hence 


f [rli^ 

Jal, J 


i^xx 4 ” ^yy'} dBody 


ri’t+i rv’l+i 

= / / (' 


^TkTk + ^BkOk)^ dTkddk- 


Now using integration by parts repeatedly we can show that 


dTkddk 


r'Pi+i rVj+i 

(2.22) / / {{<^TkTkf + ‘^ i^Tkek)^ i^Bkekf) 

Ji,f Jr,j 

f'l’i+l pj+l 2 

= / / J dTkddk + 2 ^Tkek<^ek{Vj+i,^k) 

Jipf Jrjh Jxpk 

— 2 f Urk9k^9k ddk — 2 ^ iOj^Tk^Ok dTk 

+ 2 / COrkTk^Ok irk I “^i) dTk- 

Jrjk 


ddk 


Moreover 


/ / -a; (wii + Wj/j,) dxdy 

yn* . J 

f^i+i f-Bj+i 

= - / oj{uJrkTk+<^ekek)dTkddk 

J ipk J rtj 

tp¥ Tjh Tpb 

= f f {{oJrkf + (wg^)'^) dTkddk - f ^ UJUr,('n^+i,dk) ddk 

Ji,f Jf,^ Jtl:l 

+ f ^^Tk {Vj,h) ddk- f OJUJOk {Tk, i’i+i) dTk + [ '^ ’^'^Bk {n, i’i) dTk. 

y„f 
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And this gives us the following inequality 


(2.23) 


< 


+ 


rVi+i nj+i , „ 

/ / ) dTkdOk 

J^k j^k 

i-y*. 4 4 K.+-M.)dTA 


/■A^i 

/ wWr* (^7^+1, dOk - / wwr* (r?! , 


' 'tpi ^ 'ipi 

nhl phi 

+ / ^^9^{Tk,tpi+i)dTk- OJUJB^{Tk,'lbi) dTk. 

hi 


Let tifj (Tfc, 0jt) be a set of non-conforming elements, defined on Qy the image of 
in {Tk,9k) coordinates, given by 




N N 


'^ij in, 9k) = 


n=0 m=0 


for i > 1. We shall choose u\i = 0 for all k and i. Let 

Kj] ^f+i) = {'^k, ^i+i) . 


denote the jump in u across inter element boundaries. 

We now state and prove a stability theorem which will help to motivate the stability 
theorem 2.3 on which our numerical scheme is based. 

Theorem 2.1 For the sectoral domain Cl^, the following stability estimate holds. 


(2.24) 

j=2 i=i 








“tj)] iffi+\^9k) II 


i 
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+ X! ( II I I 

i=l i=l 
N 


+ II [«i)] 


+ S (ll Mj) (^fc> ^l) ll I + II Ki) (^fc’ ^/+l) II 

J— 2 

+ E(llKA')(‘'''’'^‘)lli,w,*fH..) 


Adding a weighted combination of (2.22) and (2.23) and summing over i and j gives 

P-25) («-)E+ 2 

j=l 2=1 

< {!) + {II) + (III) + {IV) + (i/) + {VI) + {VII) + {VIII) + {IX) 


Here the terms indicated by Roman numerals given above are as follows: 


+ 1 + 


(2.26) (/) = (ul)^dnd9, 

R \ r^hi nf+i 


2K 


) ir IT i<j)e.eSdndeA , 


N-l Ik / 


w = EE ( y ^"‘ -2 (’> 5 +- 

Jh- r^i+i 

{111) = ^2/ (<„) («*„) (In ft 


i=i 


2[(uy^j-uyJ(rfc,V'f+x)dT,) 

= E fr Wj)®. Wi)r,r. 


j=l \ 

nfc 


/ Vj+I \ 

, K Rd)^r* ^4+l) j 


N-l Ik / ripf 


{VI) = -rEE 

i=l i=l 


i+1 




[«■)«•) J 
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(VII) 

{VIII) 

{IX) 


i<^) «N)r, 

i=l 

^EE i~ Ki)eJ , and 

Kj) «j)e, {n,i^l+i) dTk 

Ki) • 


i=i 

N 4-1 

Ei: 

j=l 2=1 
N 

«E 

i=i 


‘^i+i 


Now using Lemma 2.4 we can conclude that 


(2.27) 


< C 


+ 


+ 


” fn]+, rPi+i 

EE / / {<^1) ‘iniO, 

A nUi r^i+i / n2 

EEy, y 

J = 1 i=:l ''vj 

EE/. 

j = l 2=1 ^Vj 

E V'4+i) «•)" . 


where 


C = 2 m^ |max ( , (4 + 1) (V' 4+1 - V^i) 


Choosing 72 large enough and adding 



to both sides of (2.25) and then applying (2.27) and choosing K small enough we get 
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tt { ((w.). 2 («).. J + ((«y mJO 

■^4 i/ + 1* ^ 

s TttX^ 


+ 


2^ ( E IT + fT Wj)' 

Vj-1 Jv^ Jri’l 


4-1 JV 


+ EE/; [«)] (’■^.<0*/= 

i=l j=l '>vl , 

+ (//) + (///) + (IV) + (V) + (VI) + {VII) + (VIII) + {IX) . 


We shall now estimate the terms indicated by Roman numerals in the right hand side 
of (2.28) using Theorem 2.4. We begin by estimating 


l(^^)l = 


^t+1 


i‘^tj)ek {Vj+ii^k)d6k 


N-l Ik / 
j=i i=i V '’i’i 

~ fk 2 [('*^y)efc] (^ti)T*efc (^i+i>^*) 


f^i+i r 1 

■ l,k «i)rkek (vU^^k)d9k 

-'i'i 


Now by Theorem 2.4 


(2.29) 

< 

< 


rVi+i r 1 

L 2 «) (u\Sj^^,^(rt,^„e,)de, 

'>i>i 


1 [(“y)*.] 


{n 

‘y).. 

|.(^ 


2 

,VK' 

i(^?.#+l) 

2 


i.(V’?,^?+i) 


for any positive K. 
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By the trace theorem for Sobolev spaces there exists a constant M, such that 




Choosing K = 3 ^ we have 




(2.30) 


32M 

#+l 


r^t+i r 1 

/ , 2 (u*;,.),, W+1> 


< r(lnAr)H|[(uyJ(r7j;„0fc) 






And so we conclude that 


(2.31) I (7/) I < C(lnA^)2 

j=l i=l \ 

+ I [«-),.] (’>5+.. 4)'" 

N-\ Ik 




h{'PiM+\) 

+ S ^ II Wj) Il2,nj,^. • 


j=l Z =1 


Similarly we have 


iV 


(2.32) |(V^)| < C'(lnA^)^ S y|Ki)r* (^'=>'^ 1 ) 

+ I (“Dr. 






We can estimate the terms {IV ) , (V^7) , {VIII) and {IX) in the right hand side of 
(2.26) in a similar manner. Putting all these estimates together we can write (2.28) in 
the form 

j=l i=l 

r^hi rnj+i „ 

< C (In N)Uy'y2 / (^“i i) dTkdOk 

U-i i=i Jv^ 


(2.33) 
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^-1 Ik 


j=l i=l 
N h-l 






+ w: 






+ S (tl Wj) ^ 1 ) II |.(^^^^x) + II (“4 j) ^ 4 + 1 ) 11 } 

Jl. ( [“^i+i 

+ E M/ K„)(«4L(tofte*)d«i 

i=i V 


Estimating the last two terms in the same way we get the result. □ 

We shall formulate (2.33) in a more compact form. Recollect that denotes the 
circular arc with center Ak and radius p. Let B^ be the representation of this curve in 
Tk, 9k coordinates, i.e. B^ = {(ta,, Ok) \Tk = Inp, V’f < 9k < i/^^} . 

Similarly let r*, Ffc+i be the representation of Ffc and Ffc+i in (Tk,9k) coordinates. Let 
7 i be a side of Cl^j and let 7 ; be its representation in (r^, 9k) coordinates. We then have 
the following representation of (2.33). 

Consider the function . defined on C Q*'. Then the inequality 


N Ik 


(2.34) 


< 


• O • 1 ’ *lJ 

;=2 t=l 

C(lnJV)4EE||A«‘7Tt.4)l|“s;^ 


Kj=2 i=l 


fc+1 


+ E E (Il“‘llo,,. + II<1IEO 

'n=fc^.cr„n9nfc 


+ E(lli 

Ticn*’ 

Ik 


“‘iii:.,.+ii[<iiiu+ 


IIK]|lv.J 

+ EF*/ 

i=i \ •'V'? 

+ 2yJ (My^^(«y^^^^(lnp,^fe)d0fcj 
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is valid. 


The estimate (2.34) follows immediately from (2.33). 

We need to obtain a similar estimate for T*. On each j C T a function ufj (x, y) 
is defined. Now, this C and hence Qy = for some 1. We shall use fly 
and interchangeably in what follows. 

We now need to obtain an energy inequality similar to Theorem 2.1 on T^. Inte- 
grating by parts repeatedly we get 



Here s denotes arc length along dO measured from some point on it and the line 
integral is evaluated in the clockwise direction and n denotes the outward normal to 
dO, the boundary of O. Hence 



Here Bp denotes the boundary of the circle of radius p with center at A/.. Now a simple 
calculation yields 


(2.36) 


2p' 


L 


9uIn+i d 




dx 


ds 


rw+i 

= 2 / Hn+i ) (ill P. ^k) dOk 

j^f 

- 1‘" (toft h) dh 

(toft^fc) 


+ 
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Here 

(2.37) B{9,a,b) = — 2a6sin^^^ , and 

d ^ 1_^ 
dtik p dOk 


Next suppose is such that f) ^ 0. Then f) is the straight line 
joining the points D'^ and for m £ {k,k + 1} and for some 1 < j < Mm as shown 
in Fig. 2.4. Here Mm — Jm — N- Now we can show that 


(2.38) 


^ 2 f d fdu^^^\ 

^ Janrnrm % ds \ dx ) 


'anf-^^nr 

+ p^B (vr. w«).J 




^r+i 


d(7 Y) 


DV^ 


Here denotes the tangential derivative and the normal derivative along Fm. 



Figure 2.4; Points on the adjacent boundaries 
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K R f f 

2 / JL / dxdy — R / 

r 

n(-^p) 


+ Rf 

j do.’^ .. . . r\( Bis\ dn 


Similarly for with j + 1 we get 

« */,../((¥)■• (fi)’)-. 


< 


+ 


2a: 

KR 


4«/ + 

[ f (u^'?dxdy + R [ uf 

Jar*'] Jmt*' 


On 


ds. 


Let L'^ ~ {l : ^ • We can now prove the following lemma: 

Lemma 2.1 

...» i:(.'t,/((¥)’..(S)’*mi-. 
* «/,../((^)’.(T)>* 

< (/) + {II) + (///) + (IV) + (V) + {VI) + {VII) 

+ {VIII) + {IX) + {X). 

The terms indicated by Roman numerals are as follows: 

(2.42) (/) = {p^ + ^)'£ [{^’^V''{^,y)Ydxdy, 

(//) = f [ {uP+^ {x,y))^ dxdy, 

Jh. ( r'^i+i 

mi) = E|-2y_^, (<««), (top, ejds. 

- R K\JV+l) {“tV+l)rfc dOk > , 

m = p^f ^ ^ J ((“L+l)' + K^.«)^)rf^, 
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(v^) = rE 


{( E + E )/”? 

\7^canf+^nr'= 'y,cdu^+^ Dds^ 

r 1 

+ E E I > 

w=A:^^caaf+^nrm J 


dn 


ds 


(VI) = -2p^E{ 

ieL>‘ 


E^ El I [ Qy j 

[ V.canf+^ni’* -ysCdu^+^Ods^J ^ ^ ' 


dvF{^^ d fdu^+^\ 


ds 


k+l 


+ E E / «')..«“) 

’7i=*:7scanf+^nr'm 




dcJn 


h 

(VII) = 


i=l 
k+l Mm-l 


i+i 




(VIII) = - E E E K. . (»r'),„) 

m=k i=i i.aaf+'nrm.i?^?!- 

iGf 


Dni 


DV^ 


Dk ’ 

^ KA. 


(IX) = -P^’B (V-f, Ka-«).,) 

(X) = -ffB (v.f«, 


and 


jjk+i 

cf.+. 


By % we denote an arc which is a side of for I €. L'^. Here ^ denotes the radial 
derivative and ■^the tangential derivative to the circle with center at Ak and radius 
p, i.e. ^ Moreover ^ denotes the tangential derivative and ^ the normal 

derivative to the side Fjt (Fig. 2.4)- Finally Fm^i is the open subset of the straight line 
Fm between the points and and denotes ifnri = p. 

Using the estimates (2.36 - 2.40) we obtain (2.41) .□ 

Recall from (2.8) that there exists a mapping from the unit square S to 
given by 

X = Xf+'((,r,) 

y = yTHCri) 


Similarly there exists a mapping from S to 
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We now define another semi norm in terms of the transformed variables ^ and rj: 

i^.i ^ J S J 


Let 


(2-«a) = ess^sup^ (max (max (|i5»j:r‘l) .max (|D«i^'«|))) 


Then we have the following results [15] 


C 


(2.43b) |^i(^5^)lo,s — I I l^lo,nf'''^ — ^ 




(2.43c) |n(e,r7)|;,5 < ^ 




s |..|2 


J; 


Mf+' 




(2.43d) \u{^,rj)\ls < (|Mrt,oo,5l<nr + l^riU5i<nr) 

< cfH^op+i + Ko;>+0 • 


Here J^p+i denotes the Jacobian of the mapping as defined in (2.9a - 2.9b) and 




= mm 
{i,v}es 




note that we have used the bound given in (2.9b) to arrive at the above results. 
Consider the point D\ in Fig. 2.4. Then there exist two domains and on 
whose boundary lies. Let 

K+'] (Df) = <+‘ (D‘) - »r' (D‘) , 

where f) F* is traversed first if we travel along Fjt from to 
Moreover let 


[»'+*] (Gf) = (Gf) - (G‘) . 
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We now prove the following lemma. 
Lemma 2.2 


(2.44) 


\{vn) + {VIII) + {IX) + {X)\ 
32 






r)’ 


where 


Ik 


(XI) = ClnJV ^ I [(»-«) J(G*)p+ (Of) 

Vi=2 

/c + 1 ^ra — ^m.k+l / 

(2.45) + E E (([(“”■"').] (A"))'+ ([(«”«),] (BD)' 

ra—k i=2-5m,k-iri 
/c+l 

+ E p'‘ {B (i>T, (Pm) . 

m=k 

Here Pm is ifm = k and Pm is if m = k + 1. 

We first estimate one of the terms in the right hand side of (2.44). Now 


(2.46) 


p^sin^ {ipi) - (“tiv+i)rfc (“tiv-+i) } (G'f ) 


( 

(•‘WOp. (Gf) 


(Gf) 

.4 

(G‘) 

(G?) ) 


Now by Corollary 4.80 of [42] we have that if a and b are real numbers such that 
+ 6 ^ = 1 and to is a smooth function defined on such that 

w {X!*' (f,i)) (?.»))) = E E 

n=0 m—0 

then 


(2.47) 


|(am, + 6 t 0 j) (P)f < C (InN) (|w|i_„p+i + |wi 2 ,i!f+‘) ■ 
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Using (2.47) we obtain 

P^Sin^ - (^tjv+l)rfc KWOn*,} 

< (Gf))' +([«/,+.) J (G?))'} 

/| Jt |2 J_ I * -L I * 1^ I 'i 

Treating the other terms in the same way we obtain the result.Q 


We now estimate the term {IV) in (2.42). Let tu be a smooth function on Qf , 


Then 


/ w'^ds = 'Up' (p, Ok) pd9k 

= r -/>/- ( dvkde, 

Jxl^k Jp drk\v- P ) 

= r -^vPdndOk + r -2p f wwr.dndOk 

Jp u-p Jp \v-pj 


Ji>f Jp ^-p Ji,f Jp \^-PJ 

1 r'Pi+i r'Pi+i r 

< / / tt) rfcdrfcd6'fe + 2 / / \vjWrPrkdrkdek. 

P Jxp^ Jp J^’? Jp 


And so we obtain 


(2.48) 


'uPds 


'B*n9n?.iv+i 


- [—p^vu y. 


f'' 1 /'’/’i+i f'' 

/ w rkdrkd9k + - / / (w^J rkdrkd9k. 


''Pi J p 


for any a > 0. 

Hence using (2.48) we get 


{IV) =p[ (KAr+i)' + fe+i)J) 

JdnKr,,nB^^ 


(2.49) 


< f— f- ap 

V^-P 


/ (Ki^+i)L + 2 + «iv+i) Jj,) d^dy. 
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2 2 

Choose a so large that ^ ^ and choose R> ■^ + cyp+ Then combining (2.49) 

with Lemma 2.2, we have the result 


(2.50) 




leL’^ 


32 


p 


2 

2,nf+^ 


1€L>‘ 


-api J 


< ((/) + (II) + (III) + (V) + (VI)) + (XI) . 


We now obtain an estimate for the term {VI) as defined in (2.42). 
We shall estimate the first term in {VI). Now 



Let us get an upper bound on a typical element in the sum, in the right hand side of 
the above inequality, which is of the form 


(2.52) 




dy ) ds \ dx J J 


ds 


We shall assume, to be specific, that is the image of the side ^ = 1 of the square 

S under the mapping and 50^^ is the image of the side ^ = 0 of 5 under 

the mapping Recall from (2.16 - 2.17) that there exists i,j and k such that 

= Q.V for some i, j with j > N. Hence the mapping is the mapping j 
from S to Cli j as in (2.8). This representation is needed only ioi 1 < k < p, N < j < Jk 
and 1 < i < Ikj- Now Jk = N + 0{1) and Rj < I and hence there are a fixed number 
of Q,V for which this representation is needed even if we let W — >■ oo. As such we may 
assume 


(2.53) 


max 

i,j,k,j>N 




where the norm has been defined in (2.43a). Note that Cm is independent of N. We 
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shall impose further restrictions on Cm in the second part of this paper where we 
shall examine the accuracy of our numerical scheme. Here we shall only establish the 
stability of our scheme and for that an estimate of the type (2.53) is adequate. 

Now 


(2.54a) 

(2.54b) 


We have that 




= + {um^)^ Vx, and 


X 

y 




,0 <^<l, 0<r?<l. 


Let iC, v ) ) % (?) y ) ) (?i ^) and T]y (?, rf) be the unique polynomials in ^ and r} which 

are the orthogonal projections of ?* (?, 77) , rjx (?, rf ) , rf) and ijy (,^, 77) into the space 
of polynomials of degree {N — 1) in each variable separately with respect to the usual 
inner product in ((0, 1 )^) , as defined in [ 14 ]. We now define approximations to the 

Q p~l”l o pd”! 

derivatives ■ and as follows. Let 

(2-55a) ^ 

(2.55b) 


Using the approximation results in [14] we have that 


(2.56) 




1 , 00,5 


< ifmW®-"’ |lf.lL,S 


Now e, = (yr+% Moreover by (2.53) 


\Ml 


1 771 , 00,5 


^ C'm 


for all j > N] 
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and by (2.9b) 


AiP^<\ Ji,j I < Aip'^ for all ; > N. 


So it is easy to see that 


(2.57) 


'si ^I 


< CN-^ 

l,oo,S 


for all My with j > N, and N large enough. A similar result holds for ^y,T]x and rjy. 
We are now in a position to prove the following lemma. 


Lemma 2.3 Let % be contained in T^. Then 


(2.58a) 


Here 

(2.58b) 

and 

(2.58c) 


2 f \ d f \] j 
^ iy, . ds\ dx J ^ 

< C {lnNf 


V dx ) 


+ 


1/2,7. 


dy J 


2 

1/2,7., 


+ 


E(K 


p+l|2 . L p+l|2 ^ 

16 VI ^ kn^+V ■ 

1=1 


' duP+Y ' 

dx J 


dvP+'^ Y 
dy ) J 


1/2,7. 


1/2,7. 


du^'^Y , , fduY^Y . X 




1/2,(0,1) 


1/2,(0,1) 


It is easy to see that 

/o.nN f du^Y d f [dup^Y d 

i -d^Ts [-^) - 1 [-^) Ts j 

( ^ i7i) } (1 , 1?) dr) 
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+ 


^ {( (S) + («m (%)) 

(^x - fx) + iVx “ %)) } (1, ??) dr). 


Hence 

(2.60) 



by the trace theorem and an inequality for fractional Sobolev spaces we obtain below 
along with (2.43a - 2.43d). The inequality is as follows. 

Let 


N N 

771=0 n =0 


defined on S. Then 


\w\ 


1/2, S 


< C il'u;| 


o,s 


\w\ 


hS 


< CN^ llrul 


2 

0,S 


by the interpolation inequality and the inverse inequality for differentiation in [14]. 
Thus for N large enough 


< 


32 


E 

1=1 


U‘. 


p+i ^ 



( 1 , 77 ) dr) 
( 0 , 77 ) dr) 


Now 
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Clearly ( 1 , 77 ) , (^^) ( 1 , 77 ) , (^^) ( 0 , 77 ) and (0, 77 ) are polyno- 

mials in 77 of degree at most 2N. Hence by Theorem 2.4 



for any K > 0. 


Now 



Using the above estimate, the trace theorem and (2.43a - 2.43d) we get 
(2.63) 




2 2 

l/,np+^ • 

1/2, (0,1) i=i 


.P+1 1 2 


Substituting the above estimates into (2.62) and choosing K small enough we can 
conclude that 



(2.64) 
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< C(lnNf 


m 

dy 


ihv) 



+ 


+ 




32 



. p+l 2 

“>71 


+ 



■ } 

1/2, (0.1) J 


2 

1 / 2 , ( 0 , 1 ) 


Hence we get the required result. □ 
Finally 



or a similar expression. 
Now 


The form of the expressions B (77) and C {rj) do not matter except that they are analytic 
functions of rj involving Yi'^^ and their derivatives at (1, 77) • Hence we can bound 
the derivatives of B and C as in (2.54a - 2.57). Let 24 ( 77 ) be the unique polynomial 
that is the orthogonal projection of A ( 77 ) into the space of polynomials of degree N — 1 
with respect to the usual norm defined on (0, 1) . We now define 

(«r')L(i.';) = and 

= B(r,)(uX),(hr,)+C(r,)XX),iX,n)- 


It is easy to prove as we did the estimate (2.58a - 2.58c) that 


(2.65a) 


2 p" 




< c(inwf||(«r'):. 


4. — 
1/2,7,. 32 


E 

Ki=.l 


m 


p+i| 





2.3. STABILITY ESTIMATES 


44 


where 7^ C r m for some m e {k,k + l} . 


(2.65b) 


(uP) = (1,^) 

V ''<^’"11/2,7^ ^ 1/2, (0,1) 


Hence using the estimates (2.58a - 2.58c) and (2.65a - 2.65b) we obtain 


(2.66a) 


|(V/)| < (XII) + Y,jT, 1“^' 


IQLk i=i 


where 


(XII) = C(lnJVf I Y1 ( 

i7oCT*= 




xJ 1 1 1/2,75 




(2.66b) 


+ EE E i(”r')lJ 

\m=fcjei*7»canf+inrm / > 

-VE EJ¥I(^) 


Now using Lemma 2.4 and (2.9a - 2.9b) we get 


(2.67) 


E hr' (^’2/)|o,nr‘ 

leL’’ 

4EIK"‘llL. + fEE E 


iL^ JI0.7. ^ z^ ri io,7. 

7s CT>‘ \m=k leLl^ 7,canf+^ n Tm 


Shr' (x,y)|'^,+x 

ieL>‘ 


Here the constant T is independent of N. 
Ch( jose 


R > — — h exp + (T + 1) . 

v-p 
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Adding 


I 

V EE E 


'YsCTk 


to both sides of (2.50) we obtain 


^ 1 0 , 7 s 


m=k leL>= 7sC8nf+^ fl Tri 


(2.68) X) (p"|wr'(^,2/)lo,np- 


+ ~ 1“^' 2/) I (^> l2,nf+^ 

< ((/) + (//) + (HI) + (V)) + ((XI) + (XII)) 


+ 5",^ E ill"'"] 117.-^ EE E 


7scr*= 


Tn=k leL*: .y,cdnf+^ nr„ 


H0.7s 


We now have to approximate 




A<+' = af+‘ «+') + 26^' «')i + {«r'), 


Hence 


f [ dxdy = f f m+'uj+'f d^dn, 

Jnf+^ J J{o,i)xio,i) J 


where 


£f+'m = Af+‘to« + 2Bf+‘wj, + Cf+'m„ + r>f*'me + ^r'w,, 


and Af+* = etc. Let denote the unique polynomial which is the 

orthogonal projection of into the space of polynomials of degree iV — 1 in ^ and 
r) with respect to the usual inner product in (S ) . We define , Cf '^^ and 
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Ef'^^ in the same way. 

Let 

w = + 2Bf^^w^Tj + . . . 

Then it is easy to prove that for N large enough 
(2.69) (/) = J2(^p‘ + :^'jl^^J{Au^^)^dxdy 



Substituting K = p^/2R in (II) and estimating the term (V) as before along with the 
above estimates we obtain 



Here (HI) is as defined in (2.42) and (XI) is as defined in (2.45). 

We are now in a position to prove an energy inequality for the subdomain 5^ which 
we state in the following theorem. 


Theorem 2.2 Consider the subdomain S^. Then for N large enough 

(2-71) i E E ll“« + i E iwr' «. ’!)C,ar« 


j=l i=l 

< {(1) + (2) + (3) + (4)} , 
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where 


N h 


(1) = C(lnJVf 

\j=2 i=l 

+ E (IIK]II« + IIK]|IE„ + 


.k ll|2 


1/2,7, -^11^111/2,7, 


X! X (ll“*'llo,7, II^T-fctli/2,^, 

’^=fc7!crmn9n'= 


(2) = C7(lniV)M X: (IIMIIL + Ilh: 


41 + \ui]\\l 


ytQBi 


0,7, ^ IIL“'r*J ||i/2,^, -T- ||L“«JIIi/2,7, 


(3) = C'(lnA^)^ 




7,cr* 


+ E ii[“'"]iiL+ii[(«'«):]iiE+ 1 [(»'"'): 


E E 

m=k 7scar^nrTi 


^p+i I + (^iP+i)“ 

0 , 7 ^ ^ ' CTm 


(4) = ^ (-i)”'-^^ (/B 


(P-) 


'Tscasl ^ y \ / / 


By (2.43a - 2.43d) there exists a positive constant a such that 


(2.72) 


■|«r^ '7)112,5 ^ (^>7/)ll2,i 


2,fif+^ 


for all C T* and for all k. 

Then combining (2.33) and (2.70) and using (2.72) we obtain 


N h 


(2.73) 


5EEll"‘4K0‘)lll», + |ElK"'K,-;)t 

i=l t=l ' l€Lk 

< {(1) + (3) + (4)} + {XIII) + {XIV) + {XV) . 
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Here 


{XIII) = ClnivK] ![(«*=) J{Gf) + [MJ(G^) 


,1=2 

A:+l / 

+ E E (l[(“’^‘)J(A”')f+|[K'"‘)J(or) 

m=/: i=2-<Sm,fc+i ^ 



The remaining two terms are 



Once more using Theorem 2.4 we can show that 


(2.74) l(X/K)| + l(Xr)| 


< C(liiW)^ 




2 

1/2,71 


+ 



2 

1/2, 7J 





We now consider the term {D'^)f ■ There exist tw^o domains and 

such that = 7 j and is an end point of the curve 7 /. Let us assume 

that ji Pi is the image of the mapping of the boundary 77 = 1 of S' and 

7 j P is the image of the mapping Mf of the boundary 77 = 0 of S. Further let 

Dp correspond to the image of the point ^ = 0 for both these cases. 

Now 


(2.76) 

< 


I [(«'+') jCDDf 
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Moreover (^, 1) - (^, 0) j is a polynomial in ^ of degree at most 2N. 

Now by Theorem 4.79 of [42] we have that if p (s) is polynomial of degree M defined 
on [0, 1] then 

II^^IIl°o[o,i] ^ ^ (1 + InM) ||p[li/2,[o,i] • 

Hence we obtain 

3 ( ((“r‘): «, 1) - ((. 0)) ' < if to w II [(«'^>):] II . 

Now 

- ^ (I ((“•■"')< («* " f')) (“’if + 1 ((“^'), - ^•)) (“’ Ijf) • 

Using (2.57) and the Sobolev’s embedding theorem we can conclude that 

|((«rT- (0,l)f < ({’’!)llv„- 

And as before we can show 

(WN.- «'):) (o.i)f < ^ IK’«.^)ll“,. 

Choosing N large enough we obtain 

!(«"■). - (»n:) (o.i)f < I ii«r («.»)ii’,s. 


And so we can conclude that 
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In the same way, we can conclude that 

(2.77) C\nN'£ (| [K)^] (Gf)f + | [(«“), J (GJ)Q 



Substituting (2.74), (2.76) and (2.77) into (2.73) we get the required result.D 
We have now obtained an energy inequality for any of the p subdomains S’^ into which 
we had divided our original domain. Combining these estimates we can now prove the 
main theorem of this chapter which can be interpreted as a stability estimate for the 
whole domain. 

Theorem 2.3 Consider the whole domain fi. Then for N large enough there exists a 
constant C such that 


(2.78) 
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Summing the estimate (2.71) in Theorem (2.2) over k and estimating terms as 
before the result follows. □ 

2.4 Technical Results 

In this section we prove the results which we frequently refer to in Section 2.3. 


Lemma 2.4 Let w{ 6 ) be a piecewise smooth function defined for 9 € [9i,...9m+i\ 
which has discontinuities only at the points 02 , 63 , .. . 6 m- Then 




(2.79) 


d 6 


M 


+ M (t>M+l - «l) ( w" («.) + 12 (® («/) - >" (^ 7 ))" 

J=2 


Here 


w 




w ( 6 A ) — lim w ( 6 ) . 

'' ^ '' 0<ej,9-^ej ^ ' 


Define a function s (0) as follows: 


s{ 6 ) = 


w{ 6 i) 


for 61 < 6 < 62 , 


w 


(^i) + for 0fc < 0 < 6 k+i, 2<k<M 


Then w { 6 ) may be written as 


w{ 6 ) = h ( 6 ) + s (0) 


where h ( 6 ) is a continuous function which is differentiable a.e. 
Moreover h { 61 } =0. Now 

/ w^{ 6 )d 6 < 2 [ h'^{e)d 6 + s^( 6 )d 6 j 
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Clearly 


h 



Hence 

r^M+i /w/,\2 

(0) <{e- J de. 


Prom which we can conclude that 



{^M+i — 
2 



de. 


Now 

r^M+\ * o' 

/ s\e)dd < (u;(0i)f + (u;(0i))V5](u;(0t)_^i;(0-))2 

*=2 \ j =2 , 


M 


< M (^M+i - 0i) ( {w {e,)f + ^ (u; (^t) - w {e-)f ) . 

j=2 


And so we obtain the estimate. □ 


Theorem 2.4 Let aP (s) and IP (s) be polynomials of degree P on the finite interval 
[a, . Then 


(2.80) 



< C In P 




11^11* 




Here denotes the fractional Sobolev norm on (f2) as defined in [23], when 
s is not an integer. 

Now for any 0 < e < | we have 


(2.81) 


rP 

/ 

J a 


db^ (s) 
ds 


ds 


< o 


k-e,(.a,g) 


dhP 


ds 


-\+t,{a,0) 


since the space of infinitely differentiable functions with compact support in (a, ;d) is 
dense in (a, fi) for 0 < t < ^ by Theorem 1.4. 2. 4 of [23]. 
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Next using Theorem 1.4.4.6 of [23] we have the result that the differentiation operator 
is a continuous linear operator from Wg (a, / 3 ) to (a, P ) , except when t = ^, with 

norm proportional to 77 ^- 7 . Thus we can conclude that 

Q I 


(2.82) 


db^ 


ds 






Now by the interpolation inequality from [23] 


(2.83) 


1^11 sc 


And by the inverse inequality for differentiation in [14] 


(2.84) 




Once more by the interpolation inequality 




and from (2.84) we can conclude that 




This gives us the inverse inequality for fractional Sobolev norms 


(2.85) 




Using (2.83) and (2.85) we get 


ll^ni 




( 2 . 86 ) 
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Next it is easy to see that 


(2.87) 



<C 





Substituting the relations (2.82), (2.86) and (2.87) in (2.81) we get 


( 2 . 88 ) 



< ||a'’ 





Taking the minimum over positive e we get the required result. □ 



Chapter 3 


Computational Techniques for 
Dirichlet Problems 

3.1 Introduction 

In Chapter 2 we have introduced a set of local coordinate systems, which we shall 
refer to as modified polar coordinates, in a sectoral neighborhood of every corner and 
a global coordinate system elsewhere. These modified polar coordinate systems were 
first introduced by Kondratiev in [30]. We have also derived differentiability estimates 
for the solution of the Dirichlet problem in a polygonal domain in terms of these 
new systems of coordinates. We also proved a stability theorem for a non-conforming 
spectral element representation of the solution which in a sense replicates the a’priori 
estimates we obtain for Elliptic Boundary Value Problems. 

We now seek a solution to the Dirichlet problem as in [17] which minimizes a 
weighted norm of the residuals in the partial differential equation and a fractional 
Sobolev norm of the residuals in the boundary conditions and enforce continuity by 
adding a term which measures the jump in the function and its derivatives at inter- 
element boundaries in an appropriate fractional Sobolev norm to the functional being 
minimized. 
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To solve the minimization problem we need to solve the normal equations for the 
least-squares problem. To compute the residuals in the normal equations we do not need 
to compute mass and stiffness matrices as we have to in the FEM. Instead we collocate 
the partial differential equation on a finer grid of Legendre-Gauss-Lobatto points; then 
we apply the adjoint differential operator to these residuals and “ project^' these values 
back to the original grid. Such a treatment can obviously come only from an integration 
by parts procedure and hence leads to evaluation of terms at the boundaries. These 
terms can be evaluated by a collocation procedure and the other boundary terms which 
measures the jump in the function and its derivatives at inter-element boundaries in 
an appropriate Sobolev norm and a fractional Sobolev norm of the residuals in the 
boundary conditions can be easily evaluated. 

Moreover we show that we do not need to filter the coefficients of the differential and 
boundary operators or the data. Of course, the evaluation of these residuals on each 
processor requires the interchange of boundary values between neighboring processors. 
Hence the communication involved is quite small. Thus we can compute the residuals 
in the normal equations cheaply and efficiently. 

It has been shown in [18] that the methodology can be used to compute the residual 
in the p and h-p version of the FEM without having to compute any matrices. 

We solve the normal equations by the preconditioned conjugate gradient method 
using a block diagonal preconditioner and this is nearly optimal as the condition number 
of the preconditioned system is polylog arithmic in N, the number of processors and the 
number of degrees of freedom in each element. Finally we show that if our data is 
analytic then the error we commit is exponentially small in JV. 

The spectral element functions we use are non-conforming and hence there is no set 
of common boundary values to solve. Hence the method is highly efficient for solving 
problems with Dirichlet boundary conditions. 

We now briefly review the contents of this chapter. In Section 3.2 we describe 
the symmetric formulation of the numerical scheme and present in algorithmic form 
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the steps involved in computing the residuals in the normal equations. In Section 
3.3 we describe parallelization techniques and the construction of preconditioners for 
solving the normal equations by the preconditioned conjugate gradient method using a 
parallel computer. In Section 3.4 we obtain error estimates and show that the error in 
the numerical solution is exponentially small in N. Finally in Section 3.5 we provide 
computational results. 

In Chapter 4 we shall show how to solve the Poisson equation on a polygonal 
domain with mixed Neumann and Dirichlet boundary conditions. Finally in Chapter 5 
we generalize all our results and show how to solve an elliptic boundary value problem 
with analytic coefficients on a curvilinear polygon. 


3.2 The Numerical Scheme and Symmetric Formula- 
tion 

In Chapter 2 we had divided the polygonal domain Q into p sectors . . . , Qp and 

a remaining portion We had further divided each of these subdomains into still 
smaller elements jl < i < hj, 1 < i < 1 < fc < p} . It is important to note 

here that we divide these elements into p+1 subsets jl < /fc, 1 < i < ^}\<k<p 
and |1 < f < hj, N < j < In Chapter 2 we had introduced a spec- 

tral element representation on each of these subdomains. To be specific, consider 
\1 < i < Ik, We map this onto the subdomain 

h) \Vj <n< Vj+v ‘^i <0k< fpi+i } • 

We map the remaining domains |1 < i < Ik,j, N < j ^ Jk, onto the 

unit square S = {(^, ??) [O < ^ < 1, 0 < p < 1} since there is a bijective mapping MF 
from S onto We now define a non-conforming spectral element representation on 
each of these subdomains as follows: 
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Let 


dk) = 0 


ior j = 1, 1 < i < Ik, 1 < k < p, and 

Nj 

i.'Tk, Ok) = Y,Yl “m,nrr0U 

m=l n=l 

for 1 < j < N, 1 < i < Ik,j, I k < p. Here 1 < Nj < N and we shall put precise 
bounds on Nj in Section 3.4. 

We shall relabel the remaining subdomains 

{fifj 11 < i < Ik,j, N <j <Jk, l<k <p} 

as \ I < I < L] . Then we define ufj on the unit square onto which is mapped 
by as 

N N 

m=l n=l 

We are now in a position to describe our numerical scheme, which is based on the 
stability theorem 2.3. 

Let = f (^, r]) , YF (^, r;)) , for 1 < fc < p, X + 1 < j < Jfc, 1 < i < 
where (^, p) € S. Let (^, p) denote the polynomial of degree 2N — 1 in ^ and p which 
is the orthogonal projection of fF (^, p) into the space of polynomials of degree 2N — 1 
with respect to the usual inner product on H^ {S) . 

Next, let the vertex Ak = {xk, Vk) • Let (rjt, 6k) = {xk + e'^'' cos 6k,yk + sin , 
for 1 < A: < p, 1 < j < iV, 1 < i < /jfc. Here p^ < Tk < Pj+i and < 6k < 

Consider first the case j ^l. We shall let FF {rk,6k) denote the polynomial of degree 
2Nj — 1 in Tk and 6k which is the orthogonal projection of FF {Tk,6k) into the space of 
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polynomials of degree 2Nj — 1 with respect to the usual inner product on H'^ j • 
Consider now the case j = 1. Let {rk,6k) =0. 

We now consider the boundary condition « = on Ffe. Let 

(n) = ffk (xk + cosipi, yk + sin tp^) 

for rfj <rk< and i > 1. Let Tfj (r^) be the orthogonal projection of lij{Tk) into 
the space of polynomials of degree 2Nj — 1 with respect to the usual inner product on 
{Vj, Vj+i) for j > 1. Let (rk) = 0. 

Let Fa: n be the image of the mapping of S onto correspond- 
ing to ^ = 0. Let (77) = Qk (0, T}) , (0, 77)) , where 0 < 77 < 1. We shall let 

It (77) denote the polynomial of degree 2N - I which is the orthogonal projection of 
It (77) with respect to the usual inner product in H^ (0, 1) . 

Finally, we have to consider the boundary condition u = gk on Let 

^2.j ('Tfc-i) = 9k cos , Vk-i + sin 


for rij~^ < Tk-i < r]j]pl and 1 < j < N. Once more we let ij,/ (Tt-i) denote the 
polynomial of degree 2Nj — 1 which is the orthogonal projection of into the 

space of polynomials of degree 2Nj — 1 with respect to the usual inner product in 
) fo" 2 < J < W Let Tl, (Tk-i) = 0. 


Our numerical scheme may now be formulated as follows: 

Find (Tfc, Ok) } i<k<p,i<i<ik,i<j<N ’ i^hJ 0) } i<k<p,N'+i-^<Jk,i<i<ik,j } mini- 

mizes the functional 


fEEElh”.^ 

\k=l i=2 i=l 



(3.1) 
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This functional we have introduced is closely related to the right hand side of the 
stability estimate stated in Theorem 2.3. In minimizing the functional 


(^k, 0k)}, , (v^ (e, r?)} . 


we seek a solution which minimizes the sum of weighted L^ norms of the residuals 
in the partial differential equation and a fractional Sobolev norm of the residuals in 
the boundary conditions and enforce continuity by adding a term which measures the 
sum of the squares of the jumps in the function and its derivatives at inter-element 
boundaries, in an appropriate Sobolev norm. 

The above method is essentially a least-squares method and we can obtain a solution 
by using preconditioned conjugate gradient techniques for solving the normal equations. 
To be able to do so we must be able to compute the residuals in the normal equations 
inexpensively. The solution we are seeking minimizes ({vY (rfc,^)} , (^,rj)}) . 

Now 


(U + sV) = (U) + 2ey* {SU -TG) + 0 (e^) 
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for ail V, where U is a. vector assembled from the values of 






{ fe > Vm) }o<i,m<v } . } 


l<k<p 


V is a vector similarly assembled and G is assembled from the data. Here S and T 
denote matrices. Thus we have to solve SU — TG = 0; and so we have to be able to 
compute SV — TG economically during the iterative process. We now show how this 
can be done. 


We first show how to compute 


L 






where we may take 1 < A; < p, 2 < j < Jk, I < i < hj- Here 




where the coeflBcients A\p ... , E^p are polynomials of degree iV} — 1, and 
are polynomials of degree Nj and is a polynomial of degree 2Nj - 1 in each variable. 


Define the new variables 




f}j) dfdrt 
- dlj) dsdt. 


Then 
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Here 


(iy)“ w = + C^jW„ + D^jW. + BfjW, + FtjW 

where the coefficients Afp ... , are polynomials of degree Nj - 1 and gF is a poly- 
nomials of degree 2Nj — 1. 

Moreover 

Nj Nj 

(s, t) = *) ’ ^ ■^)) = S X! and 

1=0 m =0 


The differential operator is obtained from the differential operator 


with analytic coefficients by choosing A^ j so that it is the orthogonal projection of 
into the space of polynomials of degree Nj — 1 with respect to the usual inner product 
in if ^ (—1, 1) . The other coefficients are similarly defined. 

In what follows we shall drop the sub and super-scripts and denote {L^jY by L etc. 
Thus 

Lvj — Awss + Bwst + Cwtt + Dws + Ewt + Fw. 


Let L* denote the formal adjoint of the differential operator L. Then 

L^w = (Au;)„ + {Bw)^t + {CwY - (Dw)^ - {Ew)^ + Fw. 


Integrating be parts we obtain 

f f Lv (Lu — 'g) dsdt 

= / [ vL^ {Lu — ^ dsdt 

A-i.i)" J 

+ f (^{Avs + Bvt + Dv){Lu-J) -v{A {Lu-‘^)^){l,t)dt 
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— I {{Avs + Bvt + Dv) {Lu — — V (A{Lu — '^) ) 

+ f {(Cvt + Ev) iLu~^-v{C {Lu -^\-v {B {Lu - J (s, 1) ds 

J {-1,1) 

- f {Lu - p) - u {C {Lu - -v{B {Lu - J (s, -1) ds. 

J {-l,!) 

We may write this as 


dt 


/ Lv {Lu — ^ dsdt 

= f ! vD {Lu - ^ dsdt 

J 

+ f ((-Az;^ + Dv) {Lu -^ — v{B {Lu - 'g))^ - v{A {Lu - ^) ) (1, t) dt 

- [ {{Avs + Dv) {Lu-^-v{B {Lu -^),-v {A {Lu - J (-1, t) 

+ vB {Lu -^{l,t)t,-vB {Lu - ?) (-1, t)t, 

+ f {{Cvt + Ev){Lu-^-v{B{Lu-^)^-v{C{Lu-^),){s,l)ds 

- f {{Cvt + Ev) {Lu-g)-v{B {Lu -g)),-v {C {Lu - ^),) {s, -1) ds. 

Now u, V are polynomials of degree N and A, B etc. are polynomials of degree iV — 1; 
moreover ^ is a polynomial of degree 2N — 1. Hence all the integrands are polynomi- 
als of degree AN — 2 and so the integral may be exactly evaluated by the Legendre- 
Gauss-Lobatto quadrature formula with 2N + 1 points. Let tl^, . . . ,t 2 % represent the 
{2N + 1) quadrature points and Wq^,... the corresponding weights. We shall 
also denote these points by Sg^, . . . , Sgj^. Let r denote Lu — g. Then we may write 

r r 2N 2N 

I I Lvrdsdt = ^ (LV) 

Ji-hlf J i=o i=0 

2N /2N \ 

+ E “-f E (»?". ‘f ) U (1. <?")-■ (1. ‘f ) 


j-0 

2N 


,i=0 


+ J] (1, {Dr - {Br)^ - (Ar) J (l, t^f) 

j=0 
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2N 


2N 


j=o \ i=0 J 

2N 

- E ” (-1. *?") (Dr - (Br), - {Ar),) (-1. if) 

j=0 

2N / 2N \ 

+ E E (^r. ‘f ) c {sf, 1) r (,f , 1) 

2=0 \j=Q ) 

2N 

+ E (*“"■ 1) (®'- - (S’-). - (Cr),) (sf , 1) 

i=0 

2N / 2N \ 

- E E C (»f , -1) r (.f , -1) 


t=0 \j=0 

2N 


- E (®?". -1) (Dr - (Br), - (Cr),) (sf , -l) 


2=0 


+ {vBr (1, 1) — vBr (1, —1) + vBr (—1, —1) — vBr (—1, 1)) 


Here the matrix denotes the differentiation matrix. Thus 


dt 


2N 


j=0 


if Hs a polynomial of degree less than or equal to 2N. We may write this as follows. 
Recollect that f^j is a filtered representation of f^j. Dropping sub and super scripts g 
and g are the representation of / and / in the variables s and t. Let z denote Cu — g. 
Then 

2N 2N 

[ [ Lv{Lu-g)dsdt = 

>'(- 1 , 1 )^ J i=o j=0 

2N 2N 

+ (»f«,.A(l,tf)z(l.if)) 

2=0 j=Q 
2N 

+ Y, (i> (1’ ^D) 

i=o 

2N 2N 

- EE-W".*?") «C-4(-i,‘f)-(-i.*f)) 

i=0 j=0 
2N 

- £»fv{-l,tf) ((©z-(Bz).-(Az),)(-l,(f)) 

j=0 
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2N 2N 

+ EE” (4".!) z (4". 1)) 

i=0 j=0 
2N 

+ 1) 1)) 

i =0 

, 2Ar 2N 

-EE'’ (*?". *f ) (““Cc (sf . -1) z (sf- , -1)) 

1=0 j =0 
2N 

- J2 “^) - (^^)t) («r> -1)) 

i =0 

+ {v (1, 1) Bz (1, 1) - 1; (1, -1) Bz (1, -1) 

+ V (-1, -1) Bz (-1, -1) - V (-1, 1) Bz (-1, 1)) . 


Of course, in writing the above we commit an error. It can be argued as in [19] that 
this error is spectrally small. In fact if we assume that the boundary of the domain Q, 
is composed of analytic curves, the coefficients of the differential operator are analytic 
and the data is analytic then the error committed is exponentially small in N. 

Hence there is never any need to filter the coefficients of the differential and bound- 
ary operators or the data in any of our computations. 

Let U^ denote the vector 


and let 


U(l,+i1i+j+i = U (sf ,tf) hr0<i<N,0<j<N, 


C'p" +i)W+i = « (4". 4”) for 0 < i < 2N, 0 < J < 2N. 


Similarly 


y2N 

^(2N+l)i+j 




.2N 



Then we may write 


( [ Lv (Lu - dsdt = 

(-1,1)" J 
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where i? is a matrix such that RZ^^ is easily computed. 

We now show how to evaluate the boundary terms. For this we have to examine the 
norm of H^^^ (—1,1). Now 



2 

1/2, (-1,1) 




(‘ W - 1 {v)f 
- yf 


dxdy. 


Let I (t) be a polynomial of degree less than or equal to 2W— 1. Then {I (a;) — I (y)) / (x — y) 
is polynomial of degree less than or equal to 2^" - 1 in a; and y. And so we may define 


ii'i 


2 

1/2, (-1,1) 


2N 




1=0 


2N 2N 

+ E E 

j=0 


2N 2N 

wf 




2N 


+E(”'r)Ms(‘n 


i-O 


dl 


dt 


Thus there is a symmetric positive definite matrix H^^ such that 

2N 2N 

ii'iiWu) = EEUr)H.7((tf). 

i=0 j=0 


Now a typical boundary term will be of the form 


{Pus - Qut) (s, l)-h (s) 


2 

l/2,(-l,l) 


and its variation is given by 

2N 2N _ \ 

E ((^"* - O’") («?"' 1)) "S' (("“• - 0“<) (*""■ 1) - '■ (^D) • 

2=0 j=0 

So we need to examine the term 

2N 2N 

y: e ((">'> - O’") 1)) "S' - o».) (sf, 1) - A (sf )) . 

2=0 j=0 


21V / ~ X 

= E "S' (("“> - 0“.) (sf . 1) - h (sf )) . 


Let 
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Then 


2N 

2=0 

2N / 2N 


.2N 


2N 2N 


= (^f . 1) (»r. 1) - E E « (»?". 1) w". *?") 

i=0 \i=0 / i=0 j=0 

2Ar / 2Ar \ 2N 2N 

= E“(*f .1) ECp(*“ 1)^?" -EE”{*?^^?'') {Qisr,i)4^„jcf). 

j=0 \2=0 / 1=0 j=0 


Let 


= (Vus - Qut - h) (sf, 1) , 


where P and Q are filtered representations of V and Q. Similarly h is a filtered repre- 
sentation of the actual boundary data h. Then we may write 

2N 2N 

~ (S?^, 1) - ^ (sj^)) 

1=0 j =0 


2N 


' 2N 


j=o 


2N 


. 2=0 


2N 2N 

1=0 j=Q 


- /'T/2iVNi 


(v^^yTx 


2N 


Here T is a {2N + 1)^ x {2N + 1) matrix and TX^^ can be easily computed. In 
writing the above we are again committing an error and this error can be shown to be 
exponentially small in N once more. Hence there is no need to filter the coefficients of 
the boundary operators or the data. 

Adding all the terms we obtain 



Lv {Lu — ^ dsdt 


2N 2N _ N 

+ E E - «•’<) (*<"' 1) ^ (*'")) + 

1=0 j—0 


= (RZ^^ + TX^^ + ...) = (^ 2 ^)* 0 ^^ 


+ 
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where 

O^^ = RZ^^+TX^^ + ... 

is a (2N + 1)^ vector which can be easily computed. Now there exists a matrix 
such that 

y2JV ^ QNyN 


Hence 

(■y2Ny q 2N ^ ^yNy q2N^ 

In [18] it has been shown how (G^Y can be computed. Here we just describe the 
steps involved and refer the reader to [18] for further details. 

Let be the normalizing factors used in computing the Discrete Legendre Transform 
as 

„ f 1/ (* + 1/2) 

{ 2/N 

Let {^5ij}o<i<2A/,0<jf<2iV ^ follows 

n. . - n2iv 

C't.j — '^i(2N+l)+j- 

Now we perform the following set of operations. 

1. Define Oij <r- Oij . 

2. Calculate {Aij}o<i< 2 V,o<i< 2 iv Legendre transform of {Oi,i}o<K 2 iv,o<i< 2 Ar ■ De- 
fine 

A- • 4- ■ 

^ Pi Pj 

3. Calculate j Aij /pf pf , 6<i<N,0<3<N. 


if k < N, and 
if A: = A^. 
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4. Compute the inverse Legendre trnnsjorm of /x. Define 

Si,j ^ 0 < z < 0 < ;• < AT. 

5. Define a vector J of dimension {N -{-if as 

Ji(jv+i)+i = £ij lorO < i<N, d<j <N. 

Then J = (G^Y Thus we see to compute J we do not need to compute and store 
any matrices such as the mass and stiffness matrices. 

In order to compute SV — JG we need to pass the values of Uy and its derivatives 
defined on dDy , or to its neighboring processors as well as to communicate the 
corresponding values defined on neighboring processors to the processor on which ufj is 
defined. Moreover when computing the two scalars required to update the approximate 
solution and the search direction during one step of the conjugate gradient process some 
scalars have to be passed from each processor to the root processor. The two scalars 
are then computed in the root processor and pajssed to every processor. Thus we see 
communication between processors is small. 

Finally to compute SV — JG requires O (iV^) time on a parallel computer with N 
processors since to compute the two dimensional Legendre transform on a processor 
requires 0 (N^) operations. 


3.3 Parallelization Techniques 

We define the quadratic form 




ifc=l j=2 i=l 


fc=l j=N+l i=l 
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In the same way we may dehne the quadratic form 
(3.3) V" ({»?,- {"y 

= EEg +E E E 11(^1)“”?. (e.’))ll 


A;=l j—N-hl i—l 
!|2 


2 

0,5' 


+ E (iiMiil, + ii[«iii!/,„ + 

7^cnp+i ^ 

§ § § (i ^ II J ^'=) Il3/2.(^^,^^^J 

^ Il3/2.(»jt.»7^0 

P Ik / 

+ EE (II ^i^iv+i (In A ^k) - vl^ (Inp, ^ik)||o,(^f=,,t*^j 


+ (2^t\Ar+i ) e. (In p, 6^) - (v^^) (In p, 


ll/2.(V-?.^?+i) 

+ EE (iKi (’■‘>'*f)llE(,),,j».) + ii^t^ ("■‘’'*f.+i)llU(,t,,‘„) 

k — 1 j— 1 


+ E E (ll“w(°.’))llw) + l|(''yE“’’>) 

k=l j=N+l ^ 

+ E E f||’'t.(i.>!)ii + 

k=lj=N+l V 


1 / 2 , ( 0 . 1 ) 


0 ,( 0 , 1 ) 




1/2, (0.1), 


Then by Theorem 2.3 we have 

(3-4) («• (Tk, ^fc)} . .., , {vlj (e, p)} . . ,) 

< C (In Nf ({uf,. (Tk, 9k)}. . , {v^j (^, 77 )} . 

At the same time it is easy to see that there is constant K such that 
(3.5) KV" ({«?,. (ryP.)},^,*. {«« «■’')}«*) 

2 fEEEil”y(’'‘’^‘)ll2,(,j,,5«)x(*;,v,f„)+E E ElK(f.>?)ll 


\^k — 1 j — 2 i — 1 


2 

2,5 


- >v^ {n, 9k)}i jf . , {%j ■ 
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Hence we can conclude that the two quadratic forms and are spectrally equiv- 
alent and that there exists a constant K such that 


( 3 . 6 ) 


< if (In V" 



It is clear that if we use as a preconditioner, then the condition number of the 
preconditioned system is O (In N)‘^ . 

Now is a block diagonal preconditioner where each block corresponds to the 
norm of the spectral element function defined on a particular domain. Let (^, 77) 
be the spectral element function defined on the square S to which the domain Qf j is 
mapped. Then (^, 77) is determined by its values at the points 5 ?m}o<t<isr,o<m<i\r • 
Dropping sub and superscripts we order the values of v{^i,r]m) in lexicographic order 
and denote them as 7;„ for 1 < 72 < (iV + 1)^ . Now consider the bilinear form (u, v) 
induced by the norm on S, i.e. 

(u, u) = ||Mi|ff2(5) • 


Then there is a matrix A such that 

(N+lf /(7V+1)2 

S^{u,v)= 

i=l \ i=l 

The matrix A can be determined by its columns Ad where Cj is a unit vector with a 
one in its i*'* place and zero everywhere else. 

Now using integration by parts Au can be computed in O (N^) operations in ex- 
actly the same way as we have computed the residual in the normal equations. If we 
distribute the {N ■+ if columns among the Nb processors then the matrix A can be 
computed in O (iV^) time steps on a parallel computer since Nb = O (N) . Moreover the 
L-U factorization of A can be performed in time O (N^) and stored on every processor. 
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Once this has been done the action of the inverse of the matrix on different right 
hand sides can be computed in 0 {N^) time on every processor. 

We now come to the important issue of load balancing. Since the dimension of the 
spectral elements {vk, , [ul^ is {Nj + 1)^ for 2 < y < where 

CKJ ^ Nj N and a is a positive constant, hence if we were to assign each u^j (r^, dk) 
defined on onto different processors this would cause a severe imbalance in the 
loads assigned to different processors. Alternatively we could choose Nj = N for all 
j >2 defined on QN and uA jj) defined on S. By doing so we shall be able to achieve 
perfect load balancing among individual processors, but at the cost of making many of 
the processors do extra computational work which would not increcise the accuracy of 
the numerical solution substantially. 

However we should point out that the other strategy has the drawback that the 
degree of the polynomials Nj ~ aj is data dependent since a is determined by the data. 
The strategy of choosing Nj = N is thus more robust and hence is to be preferred as 
it would apply to the most general class of data. 

Finally, since we would need to perform O {N log N) iterations to obtain the solution 
to exponential accuracy and every iteration requires time O {N^), the time required to 
compute the solution would be O {N^ log N) . 


3.4 Error Estimates 

Our treatment of error estimates is very similar to the analysis in [8]. 
From the results in Section 2.2 we have 

(3.7) f [ \^Tk^ek'^i'^ki&k)\ dr^dOk 

•'A- |al<m 


for 2 < j < N. Here 0 < fik < I- 
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Now consider a quadrilateral O'- j with j > N. Let 

where the mapping 



from S to has been defined in [8]. Then the following has been proved in Lemma 
5.1 of [8]. 

Ui j is analytic on S and for some constants C and d independent of i and k and 
|q;| = m, m = 1, 2, . . . , we have 

(3.8) K,^)| <Cm!cr 

for all j > N. 

We are now in a position to prove our main theorem which we state below. 


Theorem 3.1 Let {v^.: {Tk,0k)]. ..he a polynomial in Tk and 6k of degree Nj for all i 
and k where aj < Nj < N for some positive a for j > 2. (The choice of these constants 
will become clear in what follows) 


Let {{wj 


over all j 


that 



P N Ik 

(3.9) 

EEER- 

fc=l j=l i=l 


+ E E EIK (?.>;) II V 

Jfc=l j=N+l i=l 

< Ce-^^. 
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Here {^,r)) — u • Then using the results on approximation theory in 

Section 5 of [8] we have that there exists a polynomial (j)\j (^, rf) of degree N in each 
variable separately such that 

(3.10) 

for 0 < Z < 2 and all N > s, where Cs — Hence 

(3.11) ||I7^ - ,, < C,N-^+» (CdVf 

for j > N. 

In the same way 11^^ (r^t, 9k) = u (xk, Ok) for j < N. Then there exists a polynomial 
4>i,j {'^k, ^k) of degree Nj in Xk and Ok separately such that 

(3.12) Il(c/i <C. (IV, .)-''•+* (c^''«-’)‘'"^*’(Ai<i)'-^s-2)!) 

where 

Xk = max mp (Axp^) , i| . 

We have used (3.7) to obtain the above estimate. 

Now 

(3.13) ||C1‘. (r„ S»)||li5f _ < (C (pm")''""*’ <i”-“ ("> - 2)!)' 

Hence there exists a constant C such that 

(3.14) lie/!;, in, »i)||',sj_ < (C 


for all k and i. 
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Let fiji^,ri) be the polynomial of degree (2A^ — 1) in each variable separately 
which is the orthogonal projection of f^ j (^, rj) in (S) into the space of polynomials 
of degree 2N — 1. Then 


(3.15) 
Next, let 





F (Tfc, dk) = 6^"'*=/ (rfe, 9k) for - oo < rjk < In p, < 9k < 


Since / is an analytic function we can show that there exist constants C and d such 
that 



Now let FF {Tk, 9k) denote the polynomial of degree 2Nj - 1 in and 9k separately 
that is the orthogonal projection of PF {Tk, 9k) in {{rjj, Vj+i) x {i’ii V’i+i)) i’lto the 
space of polynomials of degree {2Nj — 1) . 

Here 


FF {Tk, 9k) = F {Tk, 9k) for <Tk< TyJ+i, i’i <h< V-f+i- 


Then we have 


(r„4) - Ffj < C,, {2N, - ir^*‘ {Cpi^r^’ t,)' 


for j > 2, where is as in (3.12). 


Finally, we have to examine the error committed in approximating the boundary 
data^ = ( 5 i,... , 5 p). 
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Now we can conclude that 

(3.17) S (C . 

Let (Tk) denote the polynomial of degree {2Nj - 1) which is the orthogonal projec- 
tion of llj (Tk) in H'^ into the space of polynomials of degree {2Nj - 1) . Here 

;• > 2 . 

Then we can conclude that 


(3.18) 


< 


'Jfc Tib 

■ 1,3 h 


,3 


2 

3/2,(t7|,r)j^+l) 


Ct, i2Nj - {Xkdf^ 


2 


Let 


t(v)=s{XS''(0,rl),Y^H(l.v)) 


for 0 < 77 < 1, for m such that F* f) ^ (j). 

Let (77) denote the polynomial of degree {2N — 1) that is the orthogonal projection 
of 4 iv) in (0, 1) into the space of polynomials of degree {2N — 1) . Then we have 


(3.19) 



< Ct {2N - 1)"^‘+® (Cd‘t!)^ 

3/2.(0,l) 


for t<2N-l. 

Now consider the set of functions ("nt, ^k)}^jk > show 


as defined in (3.1), is exponentially small in N. 
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Using the estimates we have derived it is easy to see that 


PhN 

(3.20) El = Y,YU2\^^t3('^k,0k)-Ftj{rk,ek) 

fe=l i=l j=2 
p N-1 Ik 

i.k 


o.nt- 


A;=l j=l x=l 







P N h-l 

+ EEE 

A:=l 3=1 i=l 





fc=i i=2 

A:=l j=2 

A;=l 3=2 ^ 

P ^ / / \ (l-Bi.) 

+ EE^., (2^> -!)■“'■"’ M’feo 

t 1 -* o 




3 / 2 .(nt. 77)+0 


fc=l j=2 


It remains then to estimate the terms 


^ -tt. E I 

fc=l j=zN+l i=l 

(3.21) + ^ (iiwiiy,+iiiw:iiiW+llN! 


2 

0,S 

2 

1 / 2 , 7 . 


t ib=l 2=1 ^ 

l 2 

1 1 y I. \ y- N / I Ir \ /1 /i \ 1 

+ 


+ 


I i<l>i,N+l) r, P' - {(f>i,N) r, 




l/2,(^?,V'?+l) 

2 
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+ 


+ 



2 

0 ,( 0 , 1 ) 


{4CXA0’’» - {t)jn) 


1 / 2 , ( 0 , 1 ), 


Now in Section 3.2 we had let Cfj denote the differential operator with analytic coef- 
ficients such that 




where 




We can show as in [8] that there exist constants C and d such that 


< CdTrul 


for all {^,r]) e S and 1 q:| < m. A similar statement holds for all the other coefficients 
ot£‘,. 

Now 

i^ijT “ 7 ) = + E^jUrj + F^jU. 

Here Alj is the orthogonal projection of Aij in (S) into the space of polynomials 
of degree N - 1. The other coefiicients of Lfj are similarly obtained. Therefore 

(3.22) 14, i - 4,il <VCs{N- 1)-*+" (Cd^sl) . 


The same holds true for all other coefficients. 
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Now 


2 

0,S 


\WjY ' 1‘tj - ft, 


•'IJ ■'JJ 


2 

0,5 


Hence 


(3.23) 




0,5 


< Cs {{N - (Cd^slf) + Ct [{2N - {CdHlf'^ . 


Finally we show how to estimate 




1/2, 7i 


for any 7 * C The other boundary terms can be similarly handled. Then 7 * is a 
side which is common to and for some m and n. Let us assume that 7 ^ is 
the image of the side ^ = 1 of the square S under the mapping and ^ = 0 of the 
square S under the mapping Then 

IIM1Io,7, = f dv 

Jo 

< C, (iV-2*+8 (Cd^s!)^) . 


Now 




II ((«'); (a+‘)^ + («'), («*).) (1. '!) 


2 

1/2, (0,1) 


We have 


(CC‘) Jl,„) = {e‘),+ (te‘), (rC').) (1,>?) 
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(0,5) = ((c/r')^ ({;«), + iiV\ ««),) (0,5) 


Moreover 


(£^-"').(l,5)=(f^«).(0,5) 


1 / 2 , ( 0 , 1 ) ^ 11 ^^ 111 , 00 ,( 0 , 1 ) • 11^11 1 / 2 ,( 0 , 1 ) 


Now are analytic functions of ^ and r] on S and satisfy 




for (^, 17 ) € S and |q:| <m. So we can show that 


(3.24) II (( a +')^ - (I 


) <Cs(N-^^+^{Cd^s\f) 

JxJ 1 / 2 ,( 0 , 1 ) 


1 11 , 00 ,( 0 , 1 ) 


<C\nN 


using Sobolev’s embedding theorem and (2.47). Putting all these estimates together 
we can conclude that 


(3.25) 


lll('^):ill V 2 ,,. S C. In N (Cd^s\f) . 


Similarly we can show that 


(3.26) 


< Cs (iV- 2 ®+® In AT (Cd^s!)^) 



3.4. ERROR ESTIMATES 


81 


Hence we can conclude that the term E 2 , defined in (3.21) satisfies 

(3.27) E 2 < CpI (Cs (Ar-2^+8 iniv (C'd^s!)2) + c* {2N - 1)"^‘+® (Cd^t!)^) . 

We now use Sterling’s formula 


n! ~ \/2'Kne 


to simplify this expression. We choose 


/ 

o^j < Nj < PN, where 0 < a and P <1, 
Sj = ')Nj, where 0 < 7 < 1, and 
tj = 7(2JV,-1). 


Moreover s = jN and r — 7 {2N - 1 ) . Then we obtain 


< (2x7^* 

k=l j=2 



+ CpI (iV® In N (27r7iV) + {2N - 1)® (27r7 {2N - 1)) . 


Select 7 so that {Xkd'y) < 1 for all k. Then there exists a constant b>0 such that the 
estimate 

(3.28) 


holds. 


Let { {wfj (re, dt) . « (f . »))}y, J minimize t" {{oj, {n , «») } (?, t,)} J 



3.4. ERROR ESTIMATES 


82 


over all J, Thenby (3.28) 

(3.29) '''(W3(-‘.«»)}„,.{<(«,»)}yJ see-*". 


Therefore we can conclude that 


(3.30) V" ({^J, M) - < ({.,) - («.,)},.„) 


where the functional has been defined in (3.3). 
Hence using the stability theorem 2.3 we obtain 


(3.31) 


P N Ik 




fc=l j=2 i-1 
P Jk ^k,j 


fc=l j=N+l i=l 


2 

2,S 


It is easy to show that 


(3.32) 


+ 

< 


P N Ik 


EEEIK j {Tk, 0k) - Uij (rfc, 9k) I 



E E Ell4(f.’)) 


k-l j=N+l t=l 


£l'yK.»)|| 


2 

2.5 


Ce 


-bN 


And using the above estimates along with (3.14) we obtain 


(3.33) 


EEEIK {rk,0k)-wtj{rk,9k)\\lf,k, 

k=l j=l i=l 

+ E E EIK«.’i)-< k.<5 

A:=l i=iV’4-l i=l 


< Ce-»^ 
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which is the desired result (3.9) .□ 

R6 . nj , ^ks : It is easy to show that we can define a corrected version of the spectral 
element solution so that it is conforming and an exponentially accurate solution in the 
(Q) norm and this is described in the next section. 


3.5 Computational Results 



Figure 3.1: The geometric mesh 


To verify the asymptotic nature of the results we have obtained we consider Pois- 
son’s equation on a polygonal domain with Dirichlet boundary conditions. We consider 
only a sectoral domain as shown in Fig. 3.1 and show that the mesh chosen (geometric, 
with ratio .15 which leads to optimal convergence) with Dirichlet boundary data gives 
exponential convergence. 

After having obtained the non-conforming solution on all the elements by the tech- 
niques described in previous sections we can make a correction to the spectral element 
functions {ufj ??) * so that the corrected spectral element functions {ufj (^, 

are conforming and the error between the exact solution and the corrected approxima- 
tion in the (Q) norm is exponentially small in N. We do this in two steps: 
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1 . First, we make a bilinear correction (^, t])}.. ^ so that { + aJJ {^, rj ) } . . ^ 

are continuous at the vertices of the rectangles on which they are defined. For 
this we define = 0 for all i and k. If P is a vertex of Oh, for j > 2 and 
P ^ dO!l^ then we choose (P) so that [uU + (P) = u{P), where u (P) 

denotes the average of the values of u at P. If however P G we choose 
alj (P) so that (P) = 0. 

2. Next, we make a correction so that { (ulj + alj + pfj) 

are conforming. Once again, we define = 0 for all i and k. If 7 is a side of 
for j >2 and 7 C then we choose so that (uy + + Pij) (P) = 0 for 

P G 7 . Otherwise we choose so that (uy + + P^j) (P) = (u + a) (P) for 

P G 7 . Now Pij has it’s traces defined on the sides of the square S. Moreover 
the traces of P^j are polynomials on the sides of S. Let 

Ati(?.i)) = hiO, 

= itiiri) , a-nd 

Ay(l.’j) = 'k(ri)- 

I 

We then define a lifting of P^j (^, 7 ) onto <S as follows 

Pii (^1 ^ (^1 (0 (1 “ ^) + ^3 (0 (1 + ^) + 4>2 iv ) (1 + 0 + < i>i iv) (1 - 6 ) • 

We now define the corrected set of spectral element functions as 

Sfj (f . = (“y + “y- + ^yO K’ ’>'> ® 

We present the results of our numerical simulations for a sector with sectoral angle 
^ ^ radius p = 1- We choose our data so that the solution has the form of 
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the leading singular solution u = r“sin(a^) where a = ^. We divide the sector into 
three equal subsectors and choose the geometric ratio pt - .15. Let N be the number 
of spectral elements in the radial direction and the number of degrees of freedom of 
each variable in every element. Table 3.1 shows the percentage of relative error 

Table 3.1: Re lative error in perc ent against N 


N 

\\4v.r% 

2 

.7462E+01 

3 

.3709E+01 

4 

.1407E+01 

5 

.5151E+00 

6 

.1736E+00 

7 

.5720E-01 

8 

.1830E-01 

9 

.5779E-02 

10 

.1798E-02 

11 

.5547E-03 

12 

.1696E-03 


in the energy norm, defined as = i|e||£;/||u||^ where ll.y^. stands for the energy 

norm, against N lov jj, = 0.15. Fig. 3.2 shows the Log of relative error in the energy 



Figure 3.2: Log of relative error vs. N 

norm on the scale In against N and the relationship is almost linear. 

^Ne next choose an analytic solution to examine how the error depends on number 
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of iterations. Table 3.2 shows the percentage of relative error in the energy nor: 


Table 3.2: Relative e rror in percent against the number of iterations 


iterations 

Ikllsi? % 

10 

.8411E+00 

20 

.2104E+00 

30 

.4151E-01 

40 

.7066E-02 

50 

.1264E-02 

60 

.2016E-03 

70 

.6389E-05 

80 

.8497E-07 

90 

.7693E-07 

100 

.7687E-07 

112 

.7687E-07 


against the number of iterations. Fig. 3.3 shows the Log of relative error in the energy 
norm on the scale In l|e||jg^ against the number of iterations. 



Figure 3.3; Log of relative error vs. the number of iterations 




Chapter 4 


Mixed Boundary Conditions 

4.1 Introduction 

In this chapter we continue to work with a set of local coordinate systems, which 
we shall refer to as modified polar coordinates and which were first introduced by 
Kondratiev in [30], in a sectoral neighborhood of every corner and a global coordinate 
system elsewhere. We first derive differentiability estimates for the solution of the 
elliptic boundary value problem with mixed boundary conditions similar to those in 
Chapter 2. 

We next prove a stability theorem for an essentially non- conforming spectral element 
representation, viz spectral element functions which are non-conforming except that 
their values are the same at the vertices of the elements on which they are defined. 
Unfortunately the Constant in this stability theorem degrades severely compared 
to the Constant we obtain for the stability theorem 2.3 which deals with the case 
of Dirichlet boundary conditions. We now seek a solution to the elliptic boundary 
problem as in Chapter 3 which minimizes the sum of the squares of weighted norms 
of the residuals in the partial differential equation and appropriate fractional Sobolev 
norms of the residuals in the boundary conditions and enforce continuity by adding 
a term which measures the sum of the squares in the jump of the function and its 
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derivatives at inter- element boundaries in an appropriate fractional Sobolev norm to 
the functional being minimized. If we were to substitute zero data into this functional 
then it would exactly correspond to the quadratic form in the right hand side of the 
stability theorem we have obtained. 

To solve the minimization problem we need to solve the normal equations which 
arise from the least-squares problem corresponding to the functional we are minimizing. 
We have seen in Section 3.2 that the residual in the normal equations can be computed 
cheaply and efficiently without having to compute any mass and stiffness matrices. 

We now come to the aspect of parallelization. One important difference between 
this chapter and the earlier ones is that the spectral elements {u^j] we use are not 
fully non-conforming in the sense that all the function elements have the same value 
at common vertices of • The function at the corner most elements of a vertex 

Ak is taken to be a constant. We divide the vector composed of the values of the 
spectral element functions at the Legendre-Gauss-Lobatto points into two sub vectors - 
one consisting of the common values of the spectral elements at the vertices of all the 
elements and the other consisting of the remaining values. The dimension of the first 
sub vector, which can be thought of as a vector of common boundary values, is small 
compared to FEM. 

To solve the system of normal equations we first need to solve a much smaller system 
of equations corresponding to the Schur complement of the sub vector of common 
boundary values. However to compute the residual for the system of equations we 
obtain for the Schur complement of the common boundary values we need to be able 
to compute the action of the inverse of the matrix corresponding to a restricted version 
of the quadratic form in the right hand side of the stability theorem on any given vector. 
The common boundary values are the common values of vertices of 

4 and the k different values of i} which is the same constant for all i and 
for a fixed k. This quadratic form also corresponds to the functional we minimize with 
zero data, restricted to those spectral element functions which vanish at the vertices 
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of all the elements are zero for for all i and k. Now we are able 

to prove a stability theorem for precisely this quadratic form, restricted to the space 
of spectral element functions which vanish at the vertices of all the elements and the 
Constant in this stability theorem is as good as that of the stability theorem 2.3. 

Hence the parallelization techniques of Chapter 3 immediately apply and we can 
precondition the matrix corresponding to this restricted quadratic form by a block 
diagonal matrix with the same matrix occurring repeatedly as its block diagonal el- 
ements. An approximate inverse for this matrix can be computed, as in Chapter 3. 
Thus we can compute the action of the block diagonal matrix acting on a given vector 
to exponential accuracy in N, the number of elements and the number of degrees of 
freedom in each variable of spectral element functions, using O {N (log N)) iterations 
of the preconditioned conjugate gradient method. 

We can therefore solve the system of equations corresponding to the Schur comple- 
ment S of the sub vector of common boundary values provided we can obtain a good 
approximation Sa to S. Since the dimension of S is so small, i.e. S is a iVs x iVs 
matrix where Nb = 0{N), we can compute an accurate approximation to S from it’s 
definition and solve the normal equations using only O logiV) iterations of the 
preconditioned conjugate gradient method. Finally, we show that the spectral element 
representation we have computed approximates the actual solution to exponential ac- 
curacy in a norm which is the sum of the squares of the norms of functions defined 
on individual elements. It is possible to make a correction to the approximate solu- 
tion such that the corrected solution is conforming and the global error between the 
corrected and actual solution is exponentially small in N in the (fi) . 


4.2 Differentiability Estimates 

Let Q be a polygonal domain with vertices Ai, A 2 , • • ■ , Ap and corresponding sides 
Fi, r 2 , . . . , Fp where Fj joins the points Ai-i and Ai. In addition let the angle subtended 
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at Aj be ojj. 

In this chapter we shall examine the solution to the problem 
(4.1a) Au = / fox{x,y)eQ., 

with Dirichlet boundary conditions 
(4.1b) u = Qj for {x,y) € Tj, Tj C 

and Neumann boundary conditions 
(4.1c) — = gj for {x, y) e Tj, Tj C 

Here |J F^^^ = dD and we shall assume that Ft°l ^ (j). In particular, we may assume 
that Fi C Ft°^. For the case when F[°5 = <j) the solution is indeterminate up to a constant 
and so some additional condition has to be specified for a unique solution to exist. F^^^ 
will be called the Dirichlet part of the boundary and F^^^ the Neumann part of the 
boundary. 

Let V = [j |Fj C Ff°l } and M = {j |Fj C F^^-^ } . Let {xj, Oj) denote the modified 
polar coordinates centered at the vertex Aj. Let 

~ ^ (-^j) • 


Then it can be shown as in Chapter 2 that 


(4.2) 


I- (D> ^ (- - 2)0^ 


Here 0 < ydj < 1 and Pj > 1 - Tr/wy (respectively Pj > I - T^f‘2ujj if Dirichlet and 
Neumann boundary conditions are imposed on the edges Fj, F^yi, F^+i flf^i ~ )• 
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4.3 The Numerical Scheme and Stability Estimates 

We shall briefly review the introductory portion of Section 2.3 to keep this chapter as 
self contained as possible. 

We first divide into subdomains. Thus we divide Q into p subdomains S\S^,... ,3^, 
where denotes a domain which contains the vertex Ak and no other, and on each 
we define a geometric mesh as has been done in [8]. 

Let i = 1, . . . , Jfc, z = 1, . . . , /fcj} be a partition of S'^ and let G = 

ULi 

In this section we shall derive a stability estimate for a quadratic form, closely 
related to the functional we will minimize, formulated in terms of an essentially non- 
conforming spectral element representation of the approximate solution, and which is 
similar to Theorem 2.3, except that the Constant in the stability estimate degrades 
drastically as a function of N. 

Such a degradation of the Stability Constant would cause a drastic degradation 
in the performance of our parallel schemes. To overcome this problem we make the 
function representations on diflFerent elements continuous at the vertices of the elements 
only. It is in this sense that the spectral element representations are essentially non- 
conforming. The method we are going to propose for the solution of the Poisson 
equation with mixed boundary conditions can therefore be considered as a Vertex 
Based Method. 

The techniques we are going to employ for the proof of the stability theorem are 
very similar to the ones we use to prove Theorem 2.3, though there are some differences 
too. We shall therefore highlight only the differences and then simply state the stability 
theorem, leaving it for the reader to verify it’s proof. 

Recall that we have defined 


= {(2;, y) ki <rk< < h < ^f+1 } 
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(4.3) % = {(Tfc, Ok) |r?J < Tfe < < Ok < V’f+i } 

iox 1 < i < Ik, I < j < N. 

We shall represent in the form 

(4-4) {Tk, Ok) = bk, l<i< h, 


where bh is a constant. We shall represent the remaining as 
(4-5) u^j {Tk, Ok)=Yl'H (^m,nr{^dk, 1 < i < 4, 2<j<N, 

n=l m=l 

where Nj will be specified later and Nj < N for all j. To proceed further we need 
to specify the forms of Recall that corresponds to Uy for some i,j,k with 
j > N. Hence there is a mapping also denoted M^j, with j > N which maps 

the interior of the unit squares to We let 

(4.6) «!■« (xr' (f , V) . i?*' ((, »)) = E E 


Now 



2 


dTkdffk = oo 


unless bk = 0. 

We shall therefore change the weights for the corner elements only. 

We begin to see why the Constant in the stability theorem 4.1, we will prove, must 
degrade so badly when we prove Lemma 4.1, which is an analog of Lemma 2.4; 


Lemma 4.1 Let w (t) be a piecewise smooth function of t for t € ( co, Jjjv+i] which 
has discontinuities only at the points ? 72 ,?? 3 , ,'>1n ^"^d which is a constant for t G 
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(- 00 , 772 ) • Then 
(4.7) 


< 


+ 



N 


At?) ( ( 7 ?^^,) +'^(w(r);)-w ( 7 ?- ))' 


7=2 


Here a is a positive constant. 

Here At? = At?^, where At?^ = 7?j.|.i — 7?j for any 2 < j < N. Also 


w (9 ) , and 
(^ 7 ) = J^Tn w (9) . 

^ T > 6<ej,e-^ej 

Define a function s (r) as follows: 

{ w (7?^+i) for 7?jv < r < 7?iv+i, 

w (%+i) + Ef=fc+i {Vj )-^{nt)) % < ^ < »7fe+i> l<k<N-l. 

Then w (r) may be written as 


w (r) = h (r) + s (r) , 


where h (r) is a continuous function which is differentiable a.e. Moreover h (r) is 
constant for —00 < r < 7?2. Now 


^ 77 ^ 4.1 / r'HN+i f*VN+i 

(4.8) / w^T)dT< 2 ( (r) dT+ (r) dr 

Jr?2 '^’'2 


Clearly 
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since h (tIm+i) = 0- Hence 


(r) < (r - ryjv+i) f 
Jn 


TtN+l 

m 



From which we can conclude that 



Now by the Sobolev’s embedding theorem 


(m) < c 


/ 

Jr)2 \ 


(r) + 



dr. 


And 


r (r) 

7-00 « 


Hence 



Next 


/ JJ2 rvN+i 

^2 (^) ga(r-, 72 )^ 7 - + / (t) dr 

-OO *^^2 

< (Vn+i) At? + ( ^ (iV - + 1) (%+i) + E W) - ^ (^7 ))' 1 


.A:=2 

N 


j=k+l 


\ fc =2 / 
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(4.10) 


/ (r) + / s 2 (r)dr 

J-oo Jno 


CN^ {AV) (4) - w . 


Combining (4.8 - 4.10) we obtain the required result.D 

When we compare Lemma 4.1 with Lemma 2.4 we see that the constant multiplying 
the right hand side of (4.7) grows like 0 (iV^) unlike the case of (2.79) in Lemma 2.4 
where the constant is independent of N. It is this difference in behavior which causes 
the severe degradation in the Constant in the stability theorem for elliptic problems 
with mixed Neumann and Dirichlet boundary conditions which we shall prove at the 
end of this section. ' 


Next, we need to prove a result similar to a Poincare inequality in the subdomain 
which has Dirichlet boundary conditions prescribed on at least part of its bound- 
ary, viz. Recall that is an open set described in terms of {x,y) 

coordinates. Moreover is divided into subdomains by arcs 7 ^ which are the sides 
of the quadrilateral subdomains into which is divided. 

Now is divided into p subdomains PIS'* = T*, for k = 1 , . . . ,p. Each 
of these subdomains in bounded by the circular arc B* with center at A^ and radius 
p, by portions of the sides Tk and Ffc+i and lastly by 5S* which consists of a set of 
piecewise analytic arcs. We shall denote by the image of the circular arcs B* in 
{Tk,dk) coordinates. 

Recall that there is an analytic mapping from flj^- onto S which we write as 

X = Xjj {?,-)). 
y = Ytjii.n) 
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Figure 4.1: Mesh on the whole domain 
where (^, 77 ) € S. By ufj (C v) we denote the function 

and whose form has already been described. 

Further we assume that |q:| < 2 

(4.11a) \D°‘x\,\D°^y\ < Cp, and 

(4.11b) Cip'^ < Jlj<C^p\ 

where is the Jacobian of the mapping M^j as in (2.9a - 2.9b). 

We can now prove Lemma 4.2 which is stated below: 

Lemma 4.2 There exists a constant C, independent of N, such that the following 
inequality 

k=l 73 CB* l=l 

<cfe E IK1I.+ E 

\j€Vj,canp+^f\^j i=i 
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holds. 


Now 


|/i< 

< Kj(0,r7))'d77 + ^^ Ki 


for j > AT 4- 1. Here JT (<^, rj) denotes the Jacobian of the mapping ML. Summing the 
above over iV + 1 < j < Ji we obtain 




(4-13) X] 

j=N+l t=l 

^ 11^' 11^7.+ E II Milo 

\75Cri 75 CTI j=N+l i=l 




Here as in Chapter 2 


KllL = / 

j 0 


where 7 ^ = 9f2i j p] Ti, for 77 + 1 < j < Ji- 

The terms il[M^]||o^^ are similarly defined for 7* C (details are similar to that of 
Chapter 2). The terms 


are seminorms defined as 


HA-^y)\UrL / 

'' |a|=Jfc 
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Now from (4.13) we can conclude that 


(4.14) 


Ji h,: 


E IKIIL+ E Ell“.'.(-.^)IIE. 

7.CB1 -y^cmCidT^ j=N+li=l 

^ c( E II»'IIL+ E IHIIL+ E EI“l-(-.9)il«l 

\7^cri 7,CT1 j=JV+l i=l / 


In the same way we can show that 

(^■15) Ell““lll.+ E IHIE+EE 

-ysCdT^OdT^ j=N+l i=l 


+ E Ell“y(*-v)|lo,i 


c E II-IE + E 


J 2 ^2,j 


I] hi- (^’2/) I 


73 carl 09^2 


^■=^■+1 i=l 




Combining (4.14) and (4.15) we obtain 

(4.16) EEII”‘IL+ E IKI1E+ E 

k=iy,cBk 7.car2naT3 nf+'cr'y?" 




^ E IK1IL+ E 

yrs canp+i n Ti oint(T^ U r*) 


nf+^cr^Ur^ / 


Reiterating the argument outlined in (4.13 - 4.16) we immediately conclude that (4.12) 


holds. □ 

Next using Lemma 4.1 it is easy to conclude that 

(4.11) EEf£ f 

fc=i t=i V^h ^ 


+ E 


j=2 


Uij{Tk,0k)f dTkdOk 


P h / N 


< CN^iiAvfj^fl 

L fc=l i=l 


=1 i=l \i=2 




+ IIK]fe. + Ell“L('"P.«‘)llM*f,*?„) 

^=1 \7aCn*,^(75)<00 ^ ^ 
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Recall that is an open set for all k, and Arj = maxi<fc<p { A%} . Moreover 7 * is the 
image of 7 ^ in {rkJk) coordinates. Combining (4.12) and (4.17) we finally obtain the 
required estimate 

(4.18) XU] f httif + X I, / (Kj) h)f dTkdBk 

k=l i=l \ j=2 ^ 

1=1 

U=1 i=l j=2 ''“i.i i=l 


< 




+ E _E 

/c-l 75Cf2^,^(7s)<oo 


U 


.Jfcl Il2 


0 , 7 , + X X 


fc=l ;; 


7.CS^ 


+ E ll[” 


,p+il |P 

1 1 0 , 7 . 

7sCnp+i je® 7aCrjn9np+i 


+E E 


u 


u 


,p+l 


.kl l|2 


0,7ji 


0,7s 


Now if we combine (4.18), which is the new estimate we have obtained in our treatment 
of elliptic problems with mixed boundary conditions, with the remaining estimates we 
have obtained in Section 2.3 we obtain the following stability theorem. 

Theorem 4.1 Consider the domain Q. Then for N large enough the following stability 


estimate holds 


p Ik 


(4.19a) E E l“al' + E E E IK- + E ll“7‘ 

t=l A;=l j=2 i=l l=l 

< CJV‘V" {<■«.’))}«,») ■ 


Here 

(4.19b) 




= (XXXll^^t-(^'=’^fe)|lo,nt, 

Vfc=l 3=2 i=l 

+ E .E (IIMII:,rII[<]IIErIIK1I«) 

7fCn*^,/x(7f)<oo 

+ E E E + IKIIU) 

lev k=i-i j^cfi n an*,/ 4 ( 7 s)<oo 
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+ EE E 

l€^f k=l-l 7,cr, n9n^^^(7,)<oo 


u[ 


k 2 
1 / 2 . 7 . 


) 


^ (llMllo,7i llt^T'it]lll/2,7( + II[^6'*]1 Ii/2,7,) 

k=i 

ElK'W^T^r) (e,^) 

/=i 

+ E (ll[•‘lfc.+llt«lllW + ilN 

7.cnp+i ^ 


2 

0,S 


2 

1/2,75 


+ E E 

fc€i?7^canp+inr*; 

+ yr V 

fc€Ar7^canp+inr/b 


n + 
0 , 7 . 


2 

1 / 2 . 7 . 




2 

1/2, 7s 


Here is the approximate tangential derivative and is the approx- 

imate outward normal derivative to the side Tk as described in Chapter 2. Moreover 
[(u)“] and {u)y are the approximate jumps in the x and y derivatives of u across Xs 
as defined in Chapter 2. 

Following the proof of Theorem 2.3 we first show that 


p h 

EE 




P N Ik 


k=zl i=l 
P Jk 

+ E E Eil«t «.’<)lll, 

fc=l j=N+l i=l 


k=l j=2 i=l 


< 


CN^ (V" {{ul • 


And so this gives us (4.19a). □ 

Theorem 4.1 gives us a stability estimate for the Poisson equation with mixed boundary 
condition similar to Theorem 2.3 which deals with the case of Dirichlet boundary 
conditions. However, unlike Theorem 2.3, the constant multiplying the right hand side 
of (4.19a) grows rapidly like AT^, where N is the degree of the function elements in each 
variable and is also proportional to the number of elements into which Q is divided. 
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We shall use the estimate (4.19a) to formulate a numerical scheme to obtain an 
approximate solution which minimizes a functional, closely related to the right hand 
side of (4.19a), in the remaining portion of this section. However we cannot parallelize 
our method in an optimal way using (4.19a). To achieve this we shall have to make 
use of the fact that the spectral elements we have defined are continuous at vertices 
and this is discussed in Section 4.4. 

In the remaining portion of this section we shall formulate our numerical scheme. 
Let fj'j = f (Xfj- (^, rj) , YA (^, r])) for 1 < fc < p, N+1 <j<JkM<i< hj, where 
(^, 77 ) G S. Let fij (^, p) denote the polynomial of degree (2N — 1) in ^ and rj, which is 
the orthogonal projection of j {(, rj) into the space of polynomials of degree {2N — 1) 
with respect to the usual inner product on H"^ (5) . Next, let the vertex Ak = ixk,yk) ■ 
Let Flj{Tk,0k) = e^'^’^f{xk + e^’‘COsek,yk + e^’‘sm9k), lox 1 < k < p, 2 < j < N, 

I <i<Ik- Here 77 } < Tfc < and V'f < 9k < fpi+v 

Consider the case j ^ 1. We shall let PF (r*,, 9^) denote the polynomial of degree 
{2Nj - 1) in rfc and 9k, which is the orthogonal projection of PP (rj,, 9k) into the space of 
polynomials of degree {2Nj - 1) with respect to the usual inner product on H^ • 

Now consider the boundary conditions 
u = pfc on Tfc for k eV, and |^ = on Ffe for A: € W. 

Let 

f u = pfc (a:fc + cos {ip^) , yk + e^'' sin (V'f ) ) for A: € X> 

^ [ it ^ for A: € AT, 

for T]^ <rk< 77 I+ 1 , where 2< j <N. 

Let Fjfc fl / <f>, be the image of the mapping of S onto corre- 

sponding to ^ = 0 . 

Let ll (rj) = gk (0, 77 ) , Yf^^ ( 0 , 77 )) , where 0 < 77 < 1 . We shall let ( 77 ) denote 

the polynomial of degree ( 2 iV - 1) which is the orthogonal projection of I) ( 77 ) into the 
space of polynomials of degree ( 2 iV- 1) with respect to the usual inner product on 
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H^ (0, l),foriV<i<4. 

Finally, we have to consider the boundary condition u = on for k eV, 

on Ffc fl for kEAf. 

Let 

u = gk [xk-i + cos (V't'+i) : Vk-i + sin (4~\+i)) 

for k eV, 

^ cos , Vk-i + sin ) 

for A; E A', 

for rij~^ < Tk-i < where 2 <j < N. 



Now we shall examine how to approximate ^ for f = 1, 2 by polynomials of degree 
{2Nj — 1) . For this it is enough to see how to approximate (r^) for each of the cases 
k eV and k E Af. 

Here we need restrict ourselves only to 2 < j < N. Let us first consider the case kEV. 
Now by (4.2) we have 


\u (rfc, Ok) 


inrtxW, *5) < K- 



Hence by the trace theorem for Sobolev spaces we have 




Let g^ j (Tfc) be the orthogonal projection of - a* into the space of polynomials of 
degree {2Nj - 1) on H^ {r}^, rjj+i) . 

Then using the results on approximation theory in [14] we get 




< C„,(2iV,-l) 




-2m,+8 ^ {C {Xkdr^-^ 
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where 


Afc - max |l, 1 mp: (AV»f) | 

and Cs = Ce^. Choosing rrij = jNj with | {Xkdj) < 1 and using Stirling’s formula for 
approximating rrijl we obtain 


(4.20) 




for some 6 > 0 if we choose aj < Nj < N for all 2 < j < N. Here a is a positive 
constant. 

Thus we may define (rk), the polynomial of degree {2Nj - 1) which approximates 
li,j in), to be 


{rk) = Ofc + glj {rk) 


for T]j <Tk < Vj+i- Next, suppose k € M. 
Then by the trace theorem for Sobolev spaces 








Hence 




Let llj {Tk) denote the orthogonal projection of if j- (r*) into the space of polynomials of 
degree Nj on H^ {r]j, 'Hj+i) • Choosing mj = jNj with Xkd'j < 1 we obtain the estimate 


‘I (’■») - (’■*) 


< Ce-*" 


(4.22) 
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for some 6 > 0 if we choose aj < < N. The above estimate applies for Neumann 

data. 

Our numerical scheme may now be formulated as follows: 

Find 

which minimizes the functional 
(4.23) 

P N 

= EEE 

k=l j-2 i=l 

Alir.,A:l||2 , 1|L.A:11|2 

fc=l 7/Cf2^,/i(7/)<oo 
m 

+ E E E 


2 

o,hf_. 




—1^ 

^ ‘m-fc+l 


0,74 


+ 


Tk 




1 1/2,74, 


^Bk ‘m-fc+1 


+ E E E 

meM fc=m-l;j;^Cr„ ^9^^^M(74)<oo 

+ E E (iihii«+iikiiiv2,^+iik 

fc=l 

+ E I ((0'‘»r‘ 


1/2,74 

\ 

'rfcJlll/2,7, -f- 1/2,7, 


1=1 


2 

l0,S 


74COP+1 


+ E (iiK"‘iiiL+ii[K'):]iik+ll[(^"):' 

K'):. - (?): 


+ E 

E 

kev^scanp+^ClTk ' 

+ E 

E ( 


0,74 


+ 


,p+i _ 


A:6JV7,canp+inrjfe 
over all { {v^j {n, 4 ) }. ., , {vlj i?)} . . G S^. 


2 ' 

1/2.74, 
2 

1/2,74, 


f?r 


y )n 

1/2, Is J 


As in Chapter 3 it can be shown that 

P N h 


E E E II - “) i'JllU + E II K*' - “) (-. V) 

1--— 1 7— ?! i — 1 


k=l j=l i=l 

< Ce-»^ 


(4.24) 
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for some constants C and 5 as AT -> oo. In other words, the approximate solution 
converges to the actual solution at an exponential rate of accuracy in N. 

We can make a correction to the spectral element solution so that the corrected 
solution is conforming and converges to the actual solution at an exponential rate in 
the H^ (fl) norm. 


4.4 Parallelization of the Numerical Scheme 


Let U denote the values of | {u^j {n, , {wy {^, . at the Legendre-Gauss- 

Lobatto quadrature points, except that the common values at the vertices of the quadri- 

r* “1 


lateral elements are counted only once. We may then write U as U = 


Ui 

Ub 


. Here 


Ub denotes the common values of at the vertices 

of the quadrilateral elements, and Ui denotes the values of the remaining elements or- 
dered as rows and concatenated in a consistent order of elements as done in Section 


3.2. 

Let be as defined in (4.19b). Let denote 

the space of functions which are continuous at the 

vertices of the elements on which they are defined. Moreover, let {Tk,dk) — bk for 
1 <i < Ik- We shall let V denote the values of | { 

the Legendre-Gauss-Lobatto points ordered as at the beginning of this section. Further 


we may write 
(4.25) 



Then there is a symmetric, positive definite matrix A such that 
(4.26) I { Uy (Tfc ,0k )} y ^ j (^) ^) } ij^k } 


{Vi)^ (Vb)^ 


Aji Aib 
{Aib)^ Abb 


Vr 

Vb 
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When we solve the minimization problem (4.23) we have to finally solve a system of 
equations of the form 

(4.27) AU = h, 


where h can also be written as 


(4.28) 


h = 


hi 

hs 


As in Section 3.2 we may replace fF by /A etc. in the data from which we assemble 
the vector h. In doing so we commit only an exponentially small error as has been 
argued in Chapter 3. We shall therefore omit the details of how this is done and the 
interested reader is referred to Section 3.2 for a full description of the steps involved. 

Now to solve the system 

AU = h, 


we use the block L- U factorization of A, viz. 


(4.29) 


J 0 



0 


1 



0 

s 


0 I 


where the Schur complement matrix S is defined as 
(4.30) S — Agjg — AJj^Aj/ Aib- 


Consequently, solving (4.27) based on the L-U factorization given in (4.29) reduces to 
solving the system of equations 

(4.31a) SUb = Lb, where 

(4.31b) Lb — hB — A^B^nhi, 


in which the matrix S is the Schur complement matrix described in (4.30). 
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As we can see from the definition of S the feasibility of such a process depends 
on our being able to compute A/bVb, AnVi and AbbVb for any V;, Vb cheaply and 
efficiently, and this can always be done since AV can be computed inexpensively, as 
has been described in Section 3.2. 

However in addition, to this it is imperative that we should be able to construct 
effective preconditioners for the matrix An, so that the condition number of the pre- 
conditioned system is as small as possible. If we are able to do this it will be possible 
to compute {Aii)~^ V/ efficiently using the preconditioned conjugate gradient method 
for any vector Vj. Now consider the space of functions 

^0 = {{<i in, (^, 7?)},^. 4 c 


such that the values of all these functions is zero at all the vertices of the elements 




on which they are defined. In particular, this implies that wfi{Tk,dk) = 0 for all 
i and k. Let W denote the vector of the values of these functions at the Legendre- 
Gauss-Lobatto points, ordered as at the beginning of this section. Then W has the 
form 


(4.32) 




Wi 

0 


and so 

(4.33) V" i) = WjAnW,. 


We now show that for the set of functions 
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it is possible to define a quadratic form, which consists of a sum of decoupled quadratic 
forms which is spectrally equivalent to the given quadratic form 


and the constant of equivalence is at most polylogarithmic in N. We shall define one 
- such quadratic form 


in what follows, and prove the above assertions we have made for it. 

To show this we have to prove a stability theorem similar to Theorem 2.3 for the 
quadratic form 


V" { « (rt, . {wy fe, I))}y, J 


for the case when the function elements 


{ «■ (rfc, 0k )} , {<• (e, V)}. .,^) € , 


i.e. vanish at the vertices of the elements on which they are defined. 

To do so we have to be able to prove a result similar to Lemma 2.2 and this we do 
in the following lemma. 

Lemma 4.3 Let w^j (^, rj) be a polynomial of degree N in ^ and rj separately, defined 
on the unit square S = [0, 1] x [0, 1] , and which is zero at all the vertices of the square. 
Then there exists a positive constant C such that 

(4.34) ,, < C (K- tt. ’!)|L + K "'iD • 


Consider w^j (^, rj) defined on [0, 1] x {0, 1] . 
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And so we can conclude that 




Jo 


-#(«■«) 


d(- 


Integrating the above with respect to ^ we obtain 


(4.35) 




/■' 

dwi. 

L 

,-K,o) 




< 




by the trace theorem for Sobolev spaces. 
Again 


((, v) ^ iC, 0 ) + ^ ^ r?') dr}'. 


Therefore 


I (^, ?7) f < 2 \wlj (e, 0)1^ + 27] [ 

Jo 


dw^ A 


dr]. 


Integrating the above with respect to ^ and 77 we get 

j \Kj V)\^ (^dr] < 2 ^ \w^j 0)f d^ + 


dWf.; 


dr] 


d^dr]. 


Combining the above with (4.35) we obtain the required result.D 
Clearly Lemma 4.3 applies equally well to any of the function elements wf j {Tk,6k) for 
2 < i < A'', 1 < i < /fc, 1 < A: < p, although with a constant Ck which depends on k. 
Taking the supremum over the constant Ck (as given in (4.34)) we conclude that 

(4.36) \wIj (rfc, ^A:)|o^jjfc . < C {n,0k)\l^^k^ + \^tj > 
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for all function elements with 1 < k < p, i < i < 2 < j < N. Here C, of course, 

denotes a generic constant. 

Combining the estimates (4.34) and (4.36) with the remaining estimates we have 
obtained in Section 2.3 we obtain the following theorem, which in essence, is the same 
as Theorem 2.3: 


Theorem 4.2 Let {«• , {< (^, 77 )}. belong to the space of func- 

tions which are zero at the vertices of the elements on which they are defined. 
Then the following estimate holds: 


(4.37) 


< 


E E E IKj . + E IK' ^)ll2.s 

k=l j=2 i=l 1=1 


for N large enough. 


In the above wf ^ (r*,, 6k) is taken to be identically zero for 1 < A: < p and 1 < i < 

hn 

Now as in Chapter 3 we define a quadratic form 


(4.38) 


W" («■ {Tk,«t) }ijj, , ((, n) }y,t) 


kzzl j=2 t=l 


/=1 


for all in, , «■ («, e defined on 




Recall that is a rectangle of the form {‘nji‘6j+i) ^ (V'o^h-i) i ^ 
square [0, 1] x [0, 1] . 

Now using the trace theorem for Sobolev spaces we have that there exists a constant 
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K such that 


(4.39) 


< 


k=l j=2 i=l k=l j=N+l i=l 


2 

2,S 


Hence we can conclude that the quadratic form and are spectrally equivalent 
and that there exists a constant K such that 


(4.40) 




To balance the loads on all the processors, and to avoid the complication of having to 
compute a, which is data dependent, we choose Nj = N for 2 < j < N and assign 
each u^j for j > 2 on a separate processor. Here a is the constant which has to be 
chosen such that Nj ~ aj so that the error is exponentially small in N. 

We can now use the quadratic form W^, which consists of a decoupled set of 
quadratic forms for the solution on each element, as a preconditioner for V^. We can 
do this by inverting the block diagonal matrix representation of since there all the 
matrices can be replaced by a single matrix of size + 2N — 3) by {N^ + 2N — 3) 
which occurs repetitively as these blocks. We can assemble this matrix by distributing 
the computation of its columns among the Nb processors. We can then compute the 
inverse of this matrix using L-U decomposition distributed among the Nb processors 
as has been described in Chapter 3. 

Thus from (4.40) we can conclude that if we were to compute (A//)”^ Ui using 
the preconditioned conjugate gradient method then the condition number of the pre- 
conditioned matrix would be O (IniV)^ . Hence to reduce the error in the approximate 
solution by a factor of e would require O (In N) iterations. Thus to compute (An) Ui 
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to within an accuracy of O (e-*") would require 0 {N\nN) iterations of the precou- 
ditioned conjugate gradient method. 

We now return to the steps involved in solving the system of equations 

AU = h. 

As a first step we need to solve the much smaller system of equations: 

SUb — hB, 


where 


hs ~ hB ~ A^b-^u 


and the Schur complement matrix S is defined to be 

S = Abb — A^b^ii ^ib- 

Now the dimension of the vector Ub is Nb = 0{N), which is proportional to the 
number of elements into which we have divided the original domain fi. Now to solve 
(4.31a) to an accuracy of O we need to be able to compute the residual 

(4.41) Rb^SYb-Ib 

with the same accuracy and in an efficient manner. Hence we have to be able to 
compute 

SVb = AbbVb ~ A^b-^ii ^ib^^b 


for any Vb- 

As we have already seen Asb^b, ^ibVb (Ajb) I i can be computed eco- 

nomically and with communication only among neighboring processors since this holds 
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true when we compute for K Hence the bottleneck in computing R, consists in 
computing Ajj AibVb to O (e and we have seen already that this can be done using 
0 (NlnN) Iterations of the preconditioned conjugate gradient method for computing 
Af/V/ for a given vector V/. 

Thus to compute the residual Rg as defined in (4.41) to an accuracy of O 
requires O(NlnN) iterations of the preconditioned conjugate gradient method for 
computing Aj/ {AjbVb) for a given vector Vb. 

We shall now briefly examine the complexity of the solution procedure for the h-p 
version of the FEM. The set of common boundary values for the FEM consists of the 
values of the spectral element functions at the vertices and sides of the elements. In [26] 
it has been shown that we can construct an approximation Sa to the Schur complement 
S such that the condition number k of the preconditioned system satisfies 

K <C{l + {lxiNf) 

for a constant C. 

Thus to solve (4.31a) to an accuracy of 0(e“*’^) will require 0{N\ogN) itera- 
tions of the preconditioned conjugate gradient method using Sq as a preconditioner. 
Now to compute the residual in the Schur complement to an accuracy of O (e~*^) 
requires 0 (AT) iterations of the preconditioned conjugate gradient method to compute 
{^7i (^/bFb)) [42]. 

Hence we need to perform 0 {N^ log N) iterations of the preconditioned conju- 
gate gradient method for computing {Ajl) Vj-, where Vj will vary after every sequence 
of O {NlogN) steps. To solve (4.31a) to an accuracy of O (e"'’^) will require time 
0 (W® log N) on a parallel computer with N processors. 

However we now show that it is possible to solve (4.31a) with a complexity which 
is less than that we have described above. 

Recall, to begin with, that 5 is a Nb x Nb matrix, where Nb 0{N) . Hence S 
is matrix of very small dimension compared to the usual size of the Schur complement 
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matrices that arises in parallelization techniques for FEM. Let e'' be a column vector 
of dimension Nb with a 1 in the place and zeroes everywhere else. Let 

5 * = 

Then S = 

Now using (4.19a) we can show that 

(4-42) ||‘S'”^|| < CiV^. 

Moreover, it is easy to show that 

(4-43) \\Aib\\<CN\ 

Now 

= Abb^'^ - ^Jb-^ii 

Suppose we compute an approximation to 5* using 0 logiV) iterations of the 
preconditioned conjugate gradient method for computing Aj^Vj. Then using (4.42) and 

(4.43) we can show that the error between and 5*^ would be 0 . If we were 

to do this for every I .< k < Nb then the number of iterations required for computing 
Sa would be O log N ) . We would then have computed the matrix Sa, which is 
an approximate representation of the Schur complement matrix S. In fact, these two 
matrices would not be just spectrally close but we would have 

(4.44) (l- )s<Sa<{l + ) S 

for N large enough where b and c are positive constants. Hence if we were to solve 
(4.31a) using Sa as a preconditioner for S then the error after steps would be 
0 . Since each step of this process requires 0 (NiogN) iterations for computing 

(A//)~^ (A/bVb) we would have performed 0 (N^^^ logN) steps of the preconditioned 
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conjugate graciient met hod for computing Aj^ {AjbVb), where would vary after 
0{NlogN) ste{)s, (luring t his process. 

This is the same as the number of iterations required to form the approximate rep- 
resentation of the S(’hur (xunplement, Hence we can solve (4.31a) to an exponential 
accuracy in A using onh O (A log iterations of the preconditioned conjugate gra- 
dient method for computing Wi. Thus we have indeed shown that it is possible 

to solve (4.31a) with a complexity less than that which usual methods require. 
Finally, it remains to compute Uj. Now from (4.27) we have 

(4.45) AuU, = h, - A!bUb. 


Clearly we can solve the above system of equations for t/j to an accuracy of 0 (e"^^) 
jsing O {N log N) iterations of the preconditioned conjugate gradient method to com- 
DUtc (.4//) ’ V/ for a given vecttm \/. 

Thus the overall complexity of the method we have described would require 
0 log N) iterations of the preconditioned conjugate gradient method to obtain a 
lolution to the system of equations (4.27) to an exponential accuracy. Moreover the 
ime required on a parallel computer with N processors would be 0(Ar® ®logiV'). 

Our treatment of error estimates is similar to the analysis in {8] and in Chapter 3. 
Ve do not go into furtlu'r details to avoid repetition. 

t.S Computational Results 

b verify the asymptotic error estimates and estimates of computational complexity we 
onsider Poisson’s equation on a polygonal domain with mixed boundary conditions. 
Ve consider only a sectoral domain with mixed boundary data and show that our 
lethod gives exponential convergence. 

The solution we have obtained by the techniques described is conforming only at 
le vertices. We make a correction to the spectral element functions 
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as in Chapter 3 sf) tijal the corrected spectral element functions j^U 

Tt t T 0 ijjk 

OTforminR ami f lir error between the exact solution and the corrected approximation 

in the (0) norm is exponentially small in N. 

We now present t,he results of our numerical simulations for a sector with sectoral 
angle a; = f and radiiis p = 1 . We choose our data so that the solution has the form of 
the leading order singular solution u = cos (a0) where a = ^ We divide the sector 
into three eqtial subsectors and choose the geometric ratio p = .15. The boundary 
conditions are as fol low's: 


du 



Figure 4.2: The geometric mesh 


Let N be the number of spectral elements in the radial direction and the number of 
degrees of freedom of each variable in every element. Table 4.1 shows the percentage 
of relative error IjeHg^ in the energy norm, defined as HeHg;! = lie||^/ llallB where \\.\\e 
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Taljlf 4.1: Hr 


^rror in perc ent against N 


N 


2 

.7649E+01 

3 

.3720E+01 

4 

.1443E+01 

5 

.5153E-00 

6 

.1740E-00 

7 

.5723E-01 

8 

.1831E-01 

9 

5780E-02 

10 

.1799E-02 

11 

.5561E-03 

12 

.1696E-03 


error in the energy norm on t he scale in ||e||g;^ against N and the relationship is almost 
linear. 



Figure 4.3: Log of relative error vs. N 

We now provide the results for the solution of the Schur complement matrix. Table 
4.2 shows the number of iterations against the number of elements N. Fig. 4.4 shows 
the Log of number of iterations on the scale In (Iterations) against the Log of number 
of elements In {N). Here iterations denotes the minimum number of iterations of the 
preconditioned conjugate gradient method needed to solve the problem to the maximum 
accuracy achievable. To check the asymptotic estimate that Iterations = 0 InN) , 




LOG OF ITERATIONS 


















Chapter 5 


General Elliptic Problems on 
Curvilinear Polygons 


5.1 Introduction 


In this chapter we generalize all the results we have obtained in previous chapters and 
we seek a solution to an elliptic boundary value problem where the differential operator 
satisfies the Babuska - Brezzi inf-sup condition. We solve the boundary value problem 
on a curvilinear polygon whose sides are piecewise analytic and we assume the boundary 
conditions are of mixed Neumann and Dirichlet type as in [7, 8]. 

We now briefly describe the contents of this chapter. In Section 5.2 we discuss 
function spaces and obtain differentiability estimates for the solution in modified po- 
lar coordinates. In Section 5.3 we obtain a stability theorem for an essentially non- 
conforming spectral element representation, viz. spectral element functions which are 
non-conforming except at the vertices of the element on which they are defined. Finally 
in Section 5.4 we make some concluding remarks. 
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5.2. FUNCTION SPACES AND DIFFERENTIABILITY ESTIMATES i: 

5-2 Punction Spaces and Differentiability Estimates 


Let 0 be a curvilinear polygon with vertices and corresponding sides 

Fi, r2, - . • p where I » joins the points and Ai . We shall assume that the sides 
Fj are analytic arcs, i.e. 

el =:[-!, 1]} 

with (fi {^) and (^) being analytic functions on I and + WiiOf >oc>0. 

By Fi we mean the open arc, i.e. the image of I = (-1, 1) . 

Let the angle subtended at Aj be Uj. 

We shall denote the boundary dO, of fi by F. Further let F = Ft^l UFbl, Fl°l = 
UieDTj, Fbl = tJigy^Fi where P is a subset of the set {i|i = 1,... ,p} and Af = 
{i |i = 1, . . . ,p}\'D. Let x denote the vector x = {xi,X 2 ) • 

Let £ be a strongly elliptic operator 

2 2 

(5.1) £ (u) = - ^ {ar,s (x) UxXr + X! ^ 

r,5=l r=l 

where {^) = Or.s {^) . K (a;) , (V (a:) are analytic functions on fi and for any (Ci,6) ^ 
R and any a: € fi 

2 

(5.2) ^ > Fo (^1 + ^1) 

with pq > 0. Moreover let the bilinear form induced by the operator £ satisfy the 
inf-sup conditions. 

In this chapter we shall consider the boundary value problem 

(5.3) £u = / onfi, 

u = on F^°^, and 


du 


A 


= gbl on F^^^, 
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where denotes the usual conormal derivative wv.- 

, ’ ^ shall now define. Let A 

denote the 2 x 2 matrix whose entries are given by 


■^r,s {2;) — Gr,i (2;) 


for r, s = 1, 2. Let iV = (A^i, iVa) denote the outward 
Then is defined as follows 


^^rinal to the curve Fj for i € Af. 


(5.4) 



We shall assume that the given data / is analytic nn ?=r ,7 

^ and is analytic on every 

closed arc Ti and g-i®! is continuous on 

We need to state our regularity estimates in terrutj 1 1 . , , 

Of local variables which are 

. We first divide Q, 

into subdomains. Thus we divide Q, into p subdomains qi nn 1 r,,- , 

^ ,SP, where S* denotes 

a domain which contains the vertex A* and no other l o,- jo 

and on each S* we define a 

geometrical mesh. Let e‘ = {n‘^, J = 1. ■ ■ • . i = 1, . . . , 4^} he a partition of 6 ‘ 
and let © = lJ^_j 6 *. 

We now put some restrictions on 6 . Let (n.St) denote polar coordinates with 
center at At. Let n = Inr*. We choose p so that the curvilinear sector with sides 
Ffc and Ffc+i, center at Ak and radius p satisfies 


defined on a geometrical mesh imposed on Q as in Section 5 of [ 8 ] 


c y 


q!" ■ 




may be represented as 


(5.5) 


= {{x,y) e D\0 < Vk < p} _ 


The geometrical mesh we have imposed on is as shown in Fig. 51. 
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Figure 5.1: Mesh on whole curvilinear domain 


Let 1 < ^ < 4 be the side of the quadrilateral ^ ©. Then 


we assume that 


(5.6a) 




y = 


(5.6b) 


, , ^ = .... , _ 

0<ri<l, l = 2,A 


y = 

and that for some C > 1 and L>1 independent of i, j, k and 1. 


(5.7) 






ds^ 




<CLH\, t = l,2,. 


We shall place further restrictions on the geometric mesh we impose on later. 

Let {rk,ek) be polar coordinates with center at Ak- Then Q* is the open set bounded 
by the curvilinear arcs Ffc, Ffc+i and a portion of the circle rk = p- We subdivide Q* into 
curvilinear rectangles by drawing N circular arcs rk = ctj = pPk^^ J = 2, • • • , iV+1, 
where yUfe < 1 and /* - 1 analytic curves C 2 ,... , Ci^ whose exact form we shall prescribe 
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Figure 5.2: Mapping curvilinear corner onto a sector 


Moreov(!r Ikj < / for all k, j where / is a fixed constant. Let 

^'k+j = {{rk, Ok) = fj (rjfc) , 0 < ffc < p} , 

j = 0, 1 in a neighborhood Ak of 0*. Then the mapping 

(5.8) rjfc = pfc, $k = f ^ ,fcY [{<Pk - V’f) fi (Pk) - {h - i’u) /o (Pk)] , 

{Vxi ~ Vi ) 

where /| is analytic in tk for j = 0, 1, maps locally the cone 

{(pjfc, <f>k) |0 < Pfc < cr, i’l < < V’u } 

onto a set containing Q* as in Section 3 of {8]. The functions /* satisfy fy (0) — 'tpi , 
/* (0) = and ff (0) = 0 for j = 0, 1. It is easy to see that the mapping defined in 
(5.8) has two bounded derivatives in a neighborhood of the origin which contains the 
closure of the open set 


= {( Pifc ) h) |0 < PA : < P > ‘•Pi < <i>k< '>Pt } • 
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We choose the 4_i curves C 2 , ... ,(7^ as; 

C't l^fc {f'k,6k) = 

for ?; = 2, . . . , Ik. Here 

Let Then we choose so that 

(5-9) max (A^f) < A ^niin (A^f)^ 

for some constant A. We need another set of local variables {Tk,6k) in a neighborhood 
of O* where rjt = Inr*. In addition we need one final set of local variables {vkAk) in 
the cone 


{(Pife, h) |0 <pk<P, < ^* } , 

where Uk = In pit- 

Let = {(rjfc, P;t) |0 < rjfc < p} n 0. Then the image S^ in {vk, (f>k) variables of S^ is 
given by 

= {i^k, 4>k) |-oo < i^jfe < Inp, i)- < (j)’' < ipt} ■ 

Now the relationship between the variables {rk,Ok) and {i'kj<f>k) is given by {Tk,0k) = 
iykAk), viz. 

(5.10) Tjfc = Vk 

Hence it is easy to see that J^k {j^k,^k), the Jacobian of the above transformation, 
satisfies Ci < [J^k {i'k,<l>k)\ < C2 for all {vk,h) ^ fo^ 0 < p < p- 
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wo sh„..ld hooo that it is not necessary to choose the system of curves we 

have chosen ,o i.npose a geometric nresh on SJ. However it is necessary to choose the 
curve rr = „ as the boundary of n‘ and no other, as will become apparent in what 
follows. Any other additional set of analytic curves which imposes a geometrical mesh 
on S,'! would do equally well. However the set of curves we have chosen is. in some 
sense, the most natural as the image % of a curvilinear rectangle fij, for j > 2 in 

(nr, ^r) variables is given by a rectangle with straight lines for sides and for j = 1 is a 

S6nii infinito strip with strsiight linss for sides. 


We now state the differentiability estimates for the solution u of (5.3) which will be 
needed in this chapter. 

Let Utj (n*, <(>,.) = U (i/i, iPt) for I'k.ipkS % for j < Nmi at = u (At) . Now there is an 
analytic mapping I S ^ QA for j > N given by (f , ,) = ({,,), yf. ,)) . 

Here S is the unit square. Let (? . t,) = u (X^j (f , ,) , FA (^, ,)) . Then we can show 

as in previous chapters that 


(5.11a) 




fori < j < A: = 1, . . . ,p, 1 < * < 7^ and 


for N <j <J^,i<i<i^^l<k<p. 

With this we have obtained all the differentiability estimates we shall need to use in 
this chapter. 


5.3 Stability Estimates 


£m = - ^ (a^,, (a;) u^Xr + Z! ^ 

r,s=l 


U 


r=l 


Let 

(5.12) 
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be a stronfily elliptic operator which satisfies the inf-sup condition. 
Henc{‘ t her(' exist s a positive constant hq > 0 such that 

2 

Or, 5 (2:) > /^o (^1 + ^2) ) 

r,a=l 


for all 2 : € O. 

Let H = Hq (Q) w'here w € (f2) ilw Q (fi) and trace (re) |p[o] = 0. Consider 

the bilinear form B (u, v) defined on H x H as follows 


(5.13) 


B {u, v) 


L(t.- 


(x) Ux, Vxr + ^br (x) Ux,v -h CUV dx. 


r=l 


Then B {v, n) is a continuous mapping from H x H R and there exists a constant 
Cl such that 


(5.14) 


\B (u,u)| < Cl ll^illjyi(n) lbllHi{n) 


for all u,v ^ Hq {Q.) . 
Moreover 


(5.15a) 

(5.15b) 


: sup 

Iff r\-L.,ixU 


B{u,v) 


ill ... IUjII 


> C 2 > 0, and 


sup B (u, e) > 0 for every O^veH. 

u€H 


Then for every continuous linear functional F {v) defined on (Q) there exists unique 
uo e (fi) such that B (uo, v) = F (r;) for all v E (0) • 

Moreover, the a'priori estimate 


(5.16) 


1 j-F(^)l 

ll^ollHi(n) - C2 


holds. 
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Now c()nsid<'r f h<* folkwing mixed boundary value problem 


(5.17a) 

(5.17b) 

(5.17c) 


•Cu = / in n, 


JqU = 
7iu = 


^Irra and 




Here the conormal derivative j^u is defined as follows. 

Let r i C and let T and N denote the unit tangent vector and unit outward normal 

at a point P on Fj which we traverse in the clockwise direction. Let T = {Ti,T 2 f and 
N = {Nu N 2 y . Then 


(5.18a) 






r,s=l 


Nr^r.s 


du 

dxs 


iV*AVxU. 


In the same way we define the cotangential derivative 


(5.18b) 



Fi 


= yy 


,rA,,fe = TUv.«. 


With this introductory background we can now proceed to develop the results we 
shall need to prove the stability theorem, on which our numerical scheme will be based. 
Let 


l<k<p,N<j<Jk,I<i< hd} ■ 


We shall relabel the elements of and write 

Consider some as shown in Figure 6.3. Then is a curvilinear 

quadrUateral whose sides are analytic arcs and the boundary is traversed in the 

clockwise direction. 

Let 7 be a smooth curve and let and T denote the unit outward normal and 
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Figure 5.3: N and T vectors on analytic arcs of curvilinear domain 


tangent vectors to 7 at a point P on 7 . Let s be arc length measured from a point on 
the curve in the clockwise direction. Then the second fundamental form is given by 

(5.19) ^{i,v) = -^-T(v=^-NCv = f^V, 

where /c = db^ is the curvature of 7 at P. Clearly trace (95) = k. 

Now we need to use Theorem 3 . 1 . 1. 2 of [23]. Let u be a smooth vector field defined 
on where v — (ui,U 2 )*- Consider the restriction of v to the boundary - 
Now 7 i) U Qi) , where 7 i are the sides of with end points 

deleted and Qi are the vertices of We shall denote by vr the projection of v on 
the tangent vector T to except at the vertices where this cannot be defined. 

Similarly by vn we shall denote the component of v in the direction of N. Thus we 
have 


= vN, and 
Vt “ 


vT 
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Lemma 5.1 Lr.t u e . Then 



We shall say that a bounded open subset of with Lipschitz boundary F has a 
piecewise boundary if F = FolJFi, where 


1. Fo has zero measure (for the arc length measure ds) 

2. Fi is open in F and each point a: € Fj has a boundary as defined in 1.2.1.1 of 
{23j. Then Theorem 3.1. 1.2 of [23] may be stated as follows: 


Let O be a bounded open subset of with Lipschitz boundary F. Assume in addition 
that F is piecewise Then for all v G ($2) we have 


(5.21) 




dVr dVs 
dXs dXr 


f S — (vnVt) - 2vt-t^n\ ds - [ {(trace (iB))u^ + ®(r^T,yr)}cJs 
Jn i ds ds } Jvi 


To apply (5.21) we define the vector field 


V = AVxU 


where A is the matrix 


(^)r,s ■“ 
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dNj^ 


{5.22a) fOlu 

(5.22b) 

(5.22c) 

Hence (5.21) takes the form 
(5.23) 


d f du 




Y. ^rar,s ~ = (to^^) • N, and 

r,5=:l 


r,s=l 


[ \Ttufdx-'y [ 


dvj. dVi 
np+i dxs dxr 
4 


dx 


Y,( Uv,^)ds-±2( fl^) If^) & 

A, * ji hi ySTj^ds \dN)^ 

A/- (fduY fduY\^ 


Now by Lemma 3. 1.3.4 of [23] the following inequality holds for all ue IP (fi) 

d'^u 


^0 Y 

r,5=l 


^ ' d^u ? ^ 

r,5,fc,f=l 


dXrdXs 


dxgdxk dxrdxi 


a.e. in Q. 

Thus it follows that 


i4Y 

r,5s=l 


a.e. in Cl. Integrating we have 
d^u 


d^u 

^ 1 0 

d'^u das,i du 

dXrdXs 

- 2-^ Qx dx ^ ^ 

r,s=l T,s,k,l=l 

^^'^dxgdxk dxr dxi 


f4 


S/ 

r,s=l 


dXrdXs 


^ 2 ^ o Q _ r du \ \-^ \ d^u 


dx 


r,5=l 


Uf WUg 

dX c dXr 


dx + 32M 


r du 
in , 1 

r=l r,s=l 


dXrdXs 


where M is a common bound for all the norms of all the 
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(5.24) 




Ovt dvs 
dXs dXr 


dx + 


2(256)M^ ^ r\^ 
To “i \dXr 


2 

dx. 


Next 



Then combining (5.23 - 5.25) we obtain the result.D 

In a neighborhood of the vertex Ak we move to polar coordinates. We take a 
curvilinear rectangle which comprises part of the sectoral neighborhood Q*' of the 
vertex Ak and consider it’s image QY in (r*, 9k) variables as shown in Figure 5.4. 



Figure 5.4: n and t vectors on analytic arcs of curvilinear domain in {Tk,9k) coordinates 


As in Chapter 2 we write the differential operator Wl in modified polar coordinates, 


where 
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2^1 = a:* + e'^^cos0fc, and 
X2 = X 2 + e^'‘ sinOk . 


Here A* = (arfjXj) • 

We would like to obtain an estimate for 


' dx 


f rl I'm.. 

Let us define the new differential operator 
(5.26) 


/ 




dTkdOk- 




d 


dyr V^'^dy, 


E o / _ du 
^ 1 “ns] 

r,5=l 


Here yi — Tk and y 2 = Ok- 
Let O* denote the matrix 

(5.27a) 

and A* denot t' the matrix 


Qk 


cos$k -sm$k 

sin^fc cos^jfc 


r* — 


®l,l 

“1,2 

'Xk 

'Xk 

^ 2,1 

“2,2 


Then it can be easily shown that 

(5.27b) A^ = (0^)‘ AO^ 

Hence, since is an orthogonal matrix, we have that 

2 

ar,sVrVs > IM> {Vl + ■ 

r,5=l 


(5.28) 
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Moreover the foilowing relations hold 


(5.29a) 

(5.29b) 

(5.29c) 

(5.29d) 



20^,2 + 0(e"*), 

^^2 ~ 1 + 0 (e""*) , 

— 20^ 2 + (e^*’) and 

0(e^^) 


as Tjt -4 —CO. 

Next let 7 be a curve given by 


Xi = a;i(s) 
^2 = 2:2(5) 


where s is arc length along the curve 7. Then the curvature k at a point P on the curve 
is given by 

_ dxi (Tx 2 dx2 d^xi 

ds ds^ ds ds^ ’ 

Let 7 be the image of the curve in (yi, 2/2) coordinate given by 

yi = yii<^), 
y2 = y2{(r) 


where a is arc length along the curve 7. Then it is easy to verify that 

, . ds 

(5.30) 


da 


= e^K 


Now we can show that the curvature k of the curve 7 is given by 

dy2 


K = + 


da 


Hence 

(5.31) 


Jcl < |kI + 1 < 
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where K is a uniform constant for all the curves % C Q* 

We shall denote by t and n the unit tangent and outward normal vector at a point 
P on 7 , the boundary of except at it’s vertices where these are not defined. 



Now once more we use Theorem 3. 1.1. 2 of [23]. Clearly for ^' > 2 is a 
bounded open subset of with Lipschitz boundary T that is a piecewise C^. Thus 
r = (Utri Tt) U (Utri where ji are the sides of the open rectangle with the 
end points removed and Qi are it’s vertices. 

Now 




2 

dy. 


Here 



as defined in (5.26). 
Then for all w 


(5.33) 


we have 

^ '■ dwr dvjs 


r 9 /" dWrOWs, 

t{li ^ - g Z/ ^ “•> 
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Here and are the projections ot a, on the normal and tangent vectors n and t 

respectively. We define 


w = A'^VyU. 


Then 

(5.34a) 

(5.34b) 

(5.34c) 


'du\ 


E d du ^ 

E ~k du 

r,i=i 

_j(. du 


So (5.33) takes the form 

t 2 ^ 

(5.35) / dy - Y 

Jn^ . I ' 


dWr dw. 


dy 


Yi ^yr 

—2 "Y f ( —1 — ( ( —1 ^ da Y [k ( ( V ^ ^ 

j=l •'7j \ / A* \ A*/ j-i J y V^/ A*= \d'^) J 

j / du \ / du \ ( 1 (^5 ^ 

^ t A* ~ V^/A-fe U«^/A*i ’ 


Now using Lemma 3. 1.3.4 of [23] we obtain 


i4E 

r,5=l 


d^u ^ 
dyrdVs 


2 




dWr dWs 
dys dyr 


2 

H -2 E 

r,5,t,/=l 


d'^u da^ i du 
^‘"'^dysdyt dyr dyi 


And by (5.29a - 5.29d) there exists a constant M such that M is a common bound for 
the norms of all the of,. Hence 


(5.36) 


< 



d^u 

dyrdys 


2 

dy 


£4 


dWr dWj 

dys dyr 


dy + 


2(256) 


t4 




tL 

r=l ''“i.i 


du 

dyr 


2 

dy. 


Thus combining (5.33), (5.35) and (5.36) we get the result.D 
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2p‘ 


JM 


d f du\ 


xdT) ^ds vaivj 


ds 


in (5.32) whern ,, C ^ ^ ^ 

a smooth cnrve in p. < ,} , „iere p < p, and let P 

be a point on 7 such that P in polar coordinates has the representation (p,.P,) ^h 

Pk = P- 


Now 

(5.37) 


eS'^V^u = O^VyU 


where is the matrix defined in (5.27a), And 


(5-38) T = 0% and 

N = O^n. 

Hence 

(3-39a) A6^'^yu{p^ 


= t^A^Vyu(p'^ 



using (5.27a), (5.37) and (5.38). Here P is the image of the point P in {yi,y 2 ) coordi- 
nates. Similarly, we have 


(5.39b) 



Thus we can conclude that 





da, and 
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Figure 5.5: Common boundary between and 


Consider the boundary 7 common to and Then the following relations 
hold (Fig. 5.5): 



Now let 7 / C dQlj for some j <N and further suppose 7 C fj where j € V. Let 
n and t be the unit outward normal and tangent vectors, respectively, defined at every 
point of 7 . Then 


(5.42) 



Here a is arc length measured from the point G (Fig. 5.6) where 

(t*A^n (cr)] 




TA^n (cr) 
n^A’^n (cr) 
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Hence 



And so we can conclude that 




Figure 5.6: Point on a curvilinear boundary 


Next let 7 m C dQi j for some j > N be such that 7 m ^ Fj where j £V. Let N and 
T be the unit normal and tangent vectors, respectively, defined at every point of 7m. 
Then 


(5.44) 


(s) = S (^) {%) ^ 


where s is arc length measured from the point G as shown in Fig. 5.6. 
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Here 


h{s) 


T^AT- 

r‘Aiv 

mAN 


{T^ANf 
NtAN ’ 


So we obtain 



Now by (5.39b) we have that 



And moreover by (5.27a) and (5.38) 


(5.46 a) 
(5.46b) 


g(G) = 

h(G) = 


Recall that 
(5.47) 


r,5=l 

= Ttu + mu 


r=l 


where 


2 

^ hr (x) Uxr + c (x) u. 

r=l 


p 


2 


f |S!Jl«fdx<2p^ 
7nf+' 



dx + 2p^ 



dx. 


Hence 
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Hoiice using Lemma 5.1 we can conclude thpit • 

cluae that there is a constant C such that 


( 5 . 48 ) 


~ V f 


dXrdXs 



In the same way we have 


^ e f ^ (ur.s (a;) + X) + c (x) 

^ \ r,s=l r=l / 

X ~ W,S (^) + fx^r (y) (y) n 


ys=l 


'^r=l 


Here 


(5-«) = (»)%.+ 5* W” 

r=l 

and y = [yi^y^) = {rk,9k) for some 
Moreover the coefficients of satisfy 


= O (e"^*) for r = 1, 2, and 
c*= = 0 (e 2 ^*) 


^ Tjt — > — OO. 
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u 


dy<2([ ^u\y+ f 


Using Lemma 5.2 we can conclude that 




f iD^D^nt/du.d^,, 



Here Qf . = x (a*, q;|+i) and ^3 is a positive constant. 

In the same way we can rewrite (5.48) in (^, rj) coordinate as 

(5.52) f f f / 

|a|=2'^'^'^ \N<1 
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Here 5 is the unit square and is the differential operator £ written in 
coordinates. 

We now need to obtain estimates for the spectral element functions in the 
which we do in the following theorem. 


Theorem 5.1 The following estimate holds 


V h 


P N Ik 


(5.53) E E K. I + E E E ll«« llE + E 11“?^' (f. 

ktzl. i^l k-l j=2 4=1 l=l 


fcrsl k^l J=2 4=1 

I fc=l j=2 i=l 


+ 


E E (lIMIIw, + IIKlilw, + IIKIIIm.) 

*=1 ~i,cn’’ 


(ll^llo,7, 


0,7, + llo,7. 


+ EE E 

k=i-x-r,cda>‘r\LA%)<°° 

(iimiIo,7, + 


) 


*=1 7.CB* 


+ EE E 

leM' k=l~l y,Cda'‘ C{Vi 


(-) 

\dn)xk 


0,75 


^ f 


I |cur'(f, 1)1' <*?<'>) 


7,cnp+i 


+ E E 

lev 7,canp+’-nr; 


/ .. ..9 

du 

( I|w|Io,75 + 

w\ 


+E E 

0 , 75 / lelf ysCdQP+^f]T, 


du \ 

m)^ 




norm 


To prove the estimate (5.53) we shaU use (5.16). To do so we have to define a 
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corrected version of the spectral element functions so that it is conforming. 

Let be a set of spectral element func- 

tions € Si^. Here 5^ is the set of spectral element functions such that uY = h a 
constant for all is a polynomial of degree N in each variable for j > 2 and the 
spectral element functions are continuous at the vertices of the domains on which they 
are defined. Then there is a set of spectral element functions, as in Chapter 3, 


{ A-j {I'k, , {Aj,- (4, 


j>N,k 


eS^ 


such that the function ^p{xi,X 2 ) defined as 

+ Ay) (I'Jfc (3;i, 0 : 2 ) , (pk (a;i, X 2 )) if {xi,X 2 ) € forj < N 
Kj + Ay) if. {xi,X2),r]{xi,X2)) if {xi,X2) € Qljforj >N 


(piXi,X2) 


is a differentiable function of its arguments and (p & Hq (f^) . 
Moreover the estimate 

(5.64) t 


< 


fc=l i=l 


jev 


A;=l j=2 x=l 


k=l j=iV+l x=l 


fc=i-i7,crin9n*,fi(7»)<oo / 


^ (|INIIo, 7. "b ll[“‘'fc]llo,% + llt^'^*3)lo,7,) 

(llMlll7, + llWllo.7. + llt“^‘3llo,^ 


k=l 7,CB* 


+ E ( Mil. + 

7»cnp+i 


du 


dT 


+E E 

o,7j/ len 7,c9ni’+inri 


ll«1IS,7. + 


du 


0.7. , 


holds. 

We now explain the notation we have used in (5.54). Let dB denote an element of 
arc length in {i'k,4>k) coordinates. Then 
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Moreover if 7,, is given by 7, = paOP+i then 


du 


dT 


o,„ 4 \ KT &T 


i.e. 


Here qj. denotes the tangential derivative in {x\,X'^ variables, 

du . 


The other terms in the right hand side of (5.54) are similarly defined. 


Now consider the bilinear form 

B{ip,v) = / I Ut^s jx) fxsVxr +y^br jx) (px..V + O^V 1 dx 

\r,s=l TTi ) 

P N Ik L 


jz=l 1=1 


/=1 


Here 


B (<P, ^ o,,, (a;) ^ br {x) ^Xr'>J + C(pv I dx 

\r,s=l r=l / 


where A is a domain contained in fl and v € Hq (0) . 


Now 


B{^p,v)aP+i = f (y]ar,s{x)<PxsVa:r+Ylbr{x)(p^^V + C(pV 

‘ Vr,s=l r=l / 


/ £,(pvdx + 
Jnf+^ 


/ (^) 


Similarly if 1 < j < N" we have 




Moreover if j = 1 


^ = if, + isf, 


dx 
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since ip is a nmstant on QY 
Finally if j = iV + i we obtain 



^^vdx + 



vda + 


/ 



vds. 


For by (5.39b) 


P 




and ds = pda. Here P 
coordinates. Now 


is any point on the circnlar arc and P is ifs image in {n, 4 ) 


P N 


(5.55) = EEE^(..»)n,.E-(..aV. 

fc=I J = I 1=1 i 

= EEE^Wd.»)„,+E^(«r.»), 

k=i 7=1 i=i ^ 

/ p ^ i 

+ I E E E « Wd. ^n.,, + E S (Ar‘, a)„.« 

\k=l 7=1 1=1 ^ 

~ L + SIu^^\dxidx2 

k=i j=i 


/nf+^ 


S^cnS.Xoo-'’. 


+ E 

7aCnP+l 



aj 


»*+E E 


it=l 7aCfl| L 



vda 


+ EE E 


'du\ 

— 1 t;da 

/e^ i=i-i 7,cr,nan*,M(7>)<oo 





+ y'Y] Z] 

r=i 7,can?+*nri '''^‘ 

+ (eE^ +EEE^ 


A 

P N Ik 


\k=:l i=l 
L 


k=l j=2 i~l 


+ E^W^‘.a)nr‘ 


/=1 
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^''^liVdTkdek = 


i,ive'^^’’dTkdek. 


\k _ 
- 


.A: _ J -<1 ifr^orffc+iCrM, and 


otherwise 


Now c,-c(A,),a. constant, and c (i,. ij) is an analytic function of 


Hence 


Xx and X2. 


£ \ivdTkd 9 k < 2 ck( f x ( f v^e^'''‘dTkdet\ ' 


for N large enough. And so we obtain 


^XlxvdTkddk < e|Aji| ||t;(xi,a:2)||o.n^_^ 


where e is exponentially small in N. Now let 2 <j<N. Then 


I’^uljvdrkdOk < 


Finally 


(£uf+^) vdx < ||£nf+^ (xx, a:2)||o ^p+i \\v {2:1, a;2)|lo,fiP- 


P N Ik 


11 ^ ^ CN^\\v (a:i, 2;2)||f_f 


A;=l j=2 i=l 
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Hence 


(5.56) 


P N 


E E E ll» < Civr^ ||„ 


Now using the trace theorem for Sobolev 

P N Ik 

EEI 

A:=l j=2 2=1 


spaces we obtain 


EEE 


And so we can conclude that 


(5.57) 




Using the Cauchy Schwartz inequality in (5.55) and using (5.56) and (5.57) 
conclude that 


we can 


< Ar{EEE|K(’'t.<'.) 

U=1 j=2 i=l 


P h 


o,sfEEE4“‘i 

''J A;=l f=l 


k |2 


+ e(e . 



Ak 


dcr+ ^ 



dn/^kj 


da 


+ 


E/ 

trV‘ 



'rci,a:2)j dxidx2+ ^ 

7s CfiP+l •- 



ds 


+ E E y (^j/*+^EElAlt 

7s CFifl^nP+i ^ ^ k=l i=l 


P h 


k |2 


P N h 


+ 


EEE INs (’■>. «t)llk + E ll^r' (^, 

A;=l j=2 2=1 1=1 

• (ee lb(2^1.^2)||o,n|^_^ +EEEii'’(^‘.«*)iii,5f, 

I fc=l i=l k=l j=2 t=l 

+ E y)fi,nr^ + E E 53 

l=l k=l j=2 i=l ^ij 

+ E / «^*+E E 

7scnp+i /€.A/'7sCanp+inri j 
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Now o 6 Hi (n) and £ satisfies .fie inf-sup conditions (6.15a - 5.15b). Hence using 
(5.16), (5.54) and (5.57) we obtain 



^ii0,7a + 


+ EE E 

iet> ^~^-"i7aCr/n^i^^,M(75)<oo 
p 

+ EE (|IMIIo,,, + IIK)IIo,5. + IIK.JIIo,s,) 

*=1 TaCB* ^ 


)) 


+ 

+ 



(llNllo,7a + i|{^il2]|lo,T,,) 


7acnp+i 


+ E E 

lev 7aC3np+inri 


'I^ll0,7s + 


du 


dT 


0,7a > 


+ 


E E /(^)“‘*»+^EE(K.r+K0 

i6^'7acanp+inri ^ ^ p-i ;-i 


p Ik 

lE 

k=l i=l 


Here e is exponentially snaall in N. 

Using (5.54) and (5.56) once more we obtain the result.D 

We now define differential operators which are second order differential op- 

erators with polynomial coefficients in i/k and (j)k of degree (N — 1) such that these 
coefficients are exponentially close approximation to the coefficients of (Hij) as has 
been described in the discussion leading to the inequality (2.69). In the same way 
we define the differential operator (f|) to be a first order differential operator with 
polynomial coefficients in Uk and (j)k such that these coefficients are exponentially close 
approximations to the coefficients of 3® has been done in (2.55a - 2.55b). The 
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other approximations are similarly defined. 
From the above it is easy to conclude that 


(5.58) 




< C'iV4(J) 


where 


^ = {EEEii(^y““y(>'-..«ir 

L A:=l j=2 i=l 
P 


o.nt,- 


KJIIo.t, + 






+ E E (iiMiic,,.+ 

k=l 7,cn*^ 

+ EE E 

l€V k=l-l -i,<ZdO.<^P[Ti,n{%)<oo 

+ EE E 

k=i-i 7,can* nri,M(75)<oo 

+ E E (ii[»iiik + iiKiiiM,+iiwiio,,.) 

fc=l 7,CBfc 
L 


.7») 


dn)^. 


0,7* 




1=1 


o,s 


+ 


+ E E 

lev 7,canp+inri 


^ 0,7> + 


7sCfiP+l 

.dT, 


+E E 

0.7./ i€Ar7^cani>+inri 




0,7. 


Adding a weighted combination of (5.51), (5.52) and (5.58) and using the techniques 
and results of previous chapters it is not difficult to prove Theorem 5.2 which is stated 
below. 


Theorem 5.2 For N large enough the estimate 


< 


EE I-?, ir+EIKf^"*' +EII“r‘«.-))lU 

k=i i=i V i=2 ' / i=i 

CN^V^ ({<• (I'k, h) ^)}ij,k) 


(5.59) 
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other approximations are similarly defined. 
From the above it is easy to conclude that 


(5.58) 


P h 


/, 


N 


EE (l“«l +Ell«u(‘'^.«li:^g +Ell“r‘(f.'))li;, 

< CN'iX) 


where 


= (EEEII(^:y‘<('^».«l£» 


V, k=l j-2 i=l 
P 


+ > > I iii«iil2,,. + IKlK,,. + 


E E (III' 

k=l 7;, CO* 

+ EE E 

l£V A;=/~l 




.dnJxk 


o>% 


+ V I 

jLmm-J ,/ III..*/ \ 

k=:l—l jsQdQ^ f)Ti^fj,(%)<oo * 

+ E E (lIMIll^. + ilKlIlw, + IIKJIIo,,.) 

A:=1 7,CBJ 
L 


+ x; w^')“»r‘(e.’)) 


l=l 


0,5 


+ 


+ E E 

lev 73Cdnp+inr( 


: “ 0,7. + 


E (ii(“iiiS,,.+iiKJ*iio,„+iiKnil,.) 

7.CfiP+l 

duYf 


dT 


+ E E 

0,7./ l£//'r,ceflP+^f\Ti 


i^y 


0,7. 


Adding a weighted combination of (5.51), (5.52) and (5.58) and using the techniques 
and results of previous chapters it is not difficult to prove Theorem 5.2 which is stated 
below. 


Theorem 5.2 For N large enough the estimate 


(5.59) 


EE 


k=l i=l 

< CN^V^ 


V J=2 / 
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where 



^•4 Conclusion 

The numerical scheme is formulated in exactly the same way as in the earlier chapters. 
Construction of preconditioners and error estimates also proceed in just the same way. 

With this we have given a complete description and analysis of a /east-squares 
approach to solving elliptic boundary value problems in polygonal domains. 

The method is a vertex based method and the spectral element functions are con- 
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tinuous only at the vertices. As a result the Schur complement matrix has small 
dimension and an accurate inverse can be computed. Hence the numerical scheme has 
a computational complexity which is less then that of FEM. 

Moreover, the construction of a preconditioner for the Schur complement matrix 
is very simple unlike the case for FEM. In fact, for problems in 3 dimensions the 
construction of preconditioners for the Schur complement matrix becomes very complex 
for FEM [40]. 

The ideas in these chapters, though they deal with problems in two dimensions, 
generalize to 3 dimensions. We intend to study these problems in 3 dimensions both 
theoretically and computationally in future work. 
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